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ABSTRACT 


This  manual  presents  a  computer  program  for  calculating  the  dust  and 
air  temperature  environment  resulting  from  nuclear  weapons  detonations. 

The  program  can  accommodate  a  wide  variety  of  weapons  detonations  and  ground 
surfaces.  A  weapons  effects  subroutine  limits  battle  environment  input  data 
to  the  bare  minimum. 

Dust  concentrations  can  be  calculated  at  any  significant  range  from 
ground  zero  and  at  any  height  to  which  dust  is  raised  by  blast  wind  erosion 
of  soil  and  subsequent  air  transport.  The  soil  may  be  considered  either 
uniform  or  variable  in  the  path  of  the  blast  wave.  Soils  are  specified  by 
particle  size  and  terminal  velocity,  size  distribution,  and  erodable  depth. 
Land  use  is  accounted  for.  The  topography  is  assumed  to  be  flat  or  gently 
rolling.  Multiple  burst  dust  concentrations  can  be  computed  for  variously 
timed  detonations.  The  effect  of  prevailing  winds  is  introduced  after 
cessation  of  blast  winds. 

Air  temperature  is  treated  as  a  dust-related  phenomenon.  The  ground, 
heated  by  thermal  radiation  prior  to  arrival  of  the  air  blast,  is  eroded 
and  transported  by  the  blast  wind.  Its  heat  is  transferred  to  air, 
primarily  by  convection,  until  thermal  equilibrium  is  attained.  The  dust 
present  from  preceding  bursts  acts  as  a  heat  shield,  transmitting  only  a 
portion  of  thermal  radiation  from  succeeding  bursts. 

Originally  intended  for  the  NIKE-X  Power  Systems  design  studies,  the 
program  represents  the  culmination  of  four  years’  work  and  a  $365,000 
theoretical  and  experimental  study.  Although  maximum  benefit  can  be  obtained 


by  usir.g  this  manual  with  companion  documents  cited  in  the  Foreword,  enough 
theory  is  included  with  a  user's  guide  so  that  this  manual  is  virtually  self- 
contained.  The  computer  program  output  has  also  been  designed  as  input  to 
a  hardened  air  intake  and  exhaust  gas  systems  computer  program. 

I 
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FOREWORD 


An  investigation  of  the  dust  and  air  temperature  environment  resulting 
from  atmospheric  nuclear  weapons  detonations  was  performed  under  Contract 
No.  DACA  73-67-C-0018  for  the  Office  of  the  Chief  of  Engineers ,  Department 
of  the  Army,  Directorate  of  Military  Construction,  Advanced  Technology 
Branch,  Washington,  D.C.  20314.  Technical  direction  was  provided  by 
Mr.  John  J.  Healy  and  Mr.  Roger  L.  Lapp,  now  with  the  U.S.  Army 
Construction  Engineering  Research  Laboratory,  Champaign,  Illinois  61820. 

The  work  was  (tone  by  the  I IT  Research  Institute,  10  West  35th  Street, 
Chicago,  Illinois  60616.  The  project  engineer  was  Mr.  James  J.  Svatosh,  Jr. 
Other  personnel  who  contributed  materially  to  the  study  were:  A.  Anderson, 
Thomas  V.  Eichler,  and  Arne  H.  Wiedermann. 

This  manual,  while  quite  adequate  for  most  uses,  is  extracted  from  a 
lengthier  report,  recourse  to  which  will  be  of  occasional  benefit  during 
initial  trial  runs  and  in  case  of  question  on  details.  Access  to  the 
complete  report  can  be  obtained  by  contacting  the  U.S.  Army  Construction 
Engineering  Research  Laboratory  (CERL).  Environmental  output  of  the 
program  can  be  used  separately  or  as  input  to  another  computer  program 
described  in  Air  Blast  Attenuation,  Technical  Manuscript  S-l  (Revised) , 

CERL,  February  1971. 
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Chapter  1  (U) 
GENERAL  INTRODUCTION 


DUST  AND  AIR-TEMPERATURE  ENVIRONMENTAL  STUDIES 

The  dust  and  air  temperature  environment  generated  by  nuclear 
weapons  detonations  results  in  particle  and  thermal  loads  which  may  be 
detrimental  to  the  operating  characteristics  of  hardened  power  systems. 

This  is  especially  so  when  the  combustion  air  intake  and  exhaust  gas 
systems  operate  continuously,  open  to  the  atmosphere  both  during  and 
after  a  nuclear  attack. 

The  purpose  of  the  reported  studies  was  to:  develop  the  theory  for 
predicting  dust  concentrations  and  air  temperature,  experimentally 
investigate  empirical  factors  which  theoretically  describe  and  control 
the  process  (soil  erosion,  dust  transport,  and  dust-to-air  heat  transfer 
rate  coefficients),  write  a  computer  program  to  embody  the  developed 
theory  and  test  data,  evaluate  potential  NIKE-X  Ballistic  Missile 
Defense  System  power  plant  sites,  and  formulate  power  plant  design 
instructions  and  criteria. 

The  lengthier  final  technical  report  mentioned  in  the  Foreword  presents 
a  complete  documentation  of  the  studies  and  results.  Its  Chapter  2 

summarizes  the  results  of  the  studies.  It  is  intended  that  the  chapter 
will  facilitate  the  designer  by  giving  quick  answers  to  the  many  design 
questions  that  arise  concerning  dust  and  air-temperature  environments. 
Chapter  3  summarizes  the  studies  performed  under  this  contract.  It  is 
an  information  source  for  those  who  require  only  a  general  knowledge  or 
summary  of  the  studies.  Chapters  4,  5,  and  6  document  the  details  of 
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the.  studies  performed  to  date.  The  appendices  of  this  report  present 
the  tools  and  necessary  input  information  whereby  a  designer  may  compute 
detailed  dust  and  air-temperature  environments. 

This  manual  contains:  Chapters  5  and  6,  the  dust  and  air  tempera¬ 
ture  theory;  Appendix  A,  the  dust  and  air  temperature  computer  program, 
input-output  descriptions,  and  sample  problems;  and  Appendix  C,  the 
weapons  effects  theory.  The  theoretical  material  is  pic-sented  as 
necessary  background  to  the  computer  program  for  the  engineer  and  is  an 
integral  part  of  what  is  usually  called  a  computer  program  writeup. 
Program  capabilities  and  limitations  become  apparent  upon  consulting 
the  developed  theory. 

One  limitation  which  was  not  thoroughly  probed  at  the  time  of 
developing  the  theory  concerns  surface  and  near  surface  bursts  and 
deserves  mention  here.  A  recent  dust  study*  indicates  that  the  dust  and 
air  temperature  computer  program  should  be  modified  to  include  the 
effect  of  upsweep,  the  lifting  of  the  hot  fireball  gases  near  ground 
zero.  In  connection  with  this,  the  effect  of  the  blast  wave  shape, 
the  particle  size  distribution,  and  the  erosion  depth  limit  need  to  be 
examined,  along  with  the  interrelation  of  these  variables.  There  appear 
to  be  surface  burst  cases  where  the  present  dust  model  may  overstate  the 
peak  dust  concentration  in  the  low  overpressure  region  by  as  much  as 
a  factor  of  two. 

*Arne  H.  Wiedermann,  Study  of  Dust  Environment  in  the  Access  Tunnels  of 
the  NORAD  Facility,  by  IIT  Research  Institute,  Chicago,  Illinois,  for  the 
U.  S.  Army  Engineer  District,  Omaha,  Nebraska,  Contract  No.  DADC  45-71-C- 
0005,  November  1970. 
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The  dust  portion  of  the  computer  program  is  subroutined  to 
evaluate  four  different  variable  conditions.  The  first  subroutine 
•  treats  a  single  nuclear  weapon  attack  and  averages  the  soil 

characteristics  in  the  vicinity  of  the  hardened  power  system.  The 

I 

second  subroutine  of  the  program  treats  a  multiple  nuclear  weapon 
attack  and  also  treats  a  uniform  soil  distribution  in  the  vicinity 
of  the  site.  These  two  subroutines  are  useful  in  determining  the 
dust  environment  at  a  site  when  particular  details  of  the  local 
terrain  are  unknown  or  the  exact  location  of  the  site  is  unknown. 

The  third  subroutine  of  the  program  treats  a  single  nuclear  weapon 
attack  and  accounts  for  the  variability  of  local  soil  characteristics 
such  as  residential  area,  grassland,  woodland,  and  cultivated  areas. 

A  fourth  subroutine  of  the  program  treats  late- time  dust  environments, 
dust  which  has  been  transported  to  the  site  because  of  ambient  surface 
winds  after  the  blast  winds  have  subsided. 

The  air  temperature  portion  of  the  computer  program  developed 
herein  is  subroutined  to  evaluate  two  variable  conditions.  The  first 
subroutine  treats  a  uniform  ground  condition  and  the  second  treats 
a  biuniform  ground  condition.  The  subroutines  include  erosion  cutoff 

« 

factors,  single-  or  multiple-burst  attacks,  and  computation  of  air 
,  temperature  as  a  function  of  both  height  and  time. 
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The  changes  made  to  the  original  IBM  7094  FORTRAN  IV  program  to 


convert  to  CDC  6000  series  FORTRAN  are  listed  below: 


CHANGES  TO  ORIGINAL  IBM  7094  FORTRAN  IV  PROGRAM 


1.  Page  A-49:  Statement  2,  formerly  for  IBM  7094  read  - 

2  FORMAT  (18A4,2I2,L2) 

This  has  been  changed  for  CDC  6600  to  read  - 
2  FORMAT  (18A1,2I2,L2) 

2.  Page  A-88:  Lines  27-31,  formerly  for  IBM  7094  read  - 

TS(N)  =  TS(N)+((C5+A(I)*T(I))*C6/EXP(E2)-(C5+A(I)*T(I1))*C7/ 

1EXP(E3))/1.5+C1*(A(I)*(3.*T(21)+C2)+B3(I))*(ERF(C4/C6)-ERF(C4/C7)) 
GO  to  106 

105  TS(N)  =  TS (N)4CGNS3*(C5+A(I)*T(I))/EXP(E1)+C1*(A(I)*(3.*T(21)+C2)+ 
1B3(I) ) *(ERF(EF1) -1 . ) 


This  has  been  changed  for  CDC  6600  to  read  - 

ERX  =  ABS(C4/C6) 

CALL  ERF(ERX»0. »ERX1»ERXA»ERXB) 

ERX  =  ABS (C4/C7) 

CALL  ERF(ERX»0. »ERX2»ERXA»ERXB) 

TS (N)  =  TS(N)+((C5+A(I)*T(I))*C6/EXP(E2)-(C5+A(I)*T(I1))*C7/ 
1EXP(E3))/1.5+C1*(A(I)*(3.*T(21)+C2)+B3(I))*(ERXL-ERX2) 

GO  to  106 
105  ERX  =  EF1 

CALL  ERF(ERX’0.  »ERX3’ERXA»ERXB) 

TS  (N)  =  TS(N)+C0NS3*(C5+A(I)*T(I))/EXP(E1)+C1*(A(I)*(3.*T(21)+C2)+ 
1B3(I))*(ERX3-1.) 


EXTRACTED 
Chapter  5  (U) 

DUST  STUDIES 

5.1  Introduction 

The  adequate  design  of  hardened  air  breathing 
facilities  requires  that  the  local  ambient  air  environment 
be  established  during  a  nuclear  attack  so  that  suitable 
compensatory  action  can  be  taken,  thus  ensuring  proper 
functioning  of  the  power  system.  The  air  in  the  vicinity 
of  the  facility  may  become  dust  ladened,  and  will  experience 
a  temperature  and  pressure  transient. 

The  majority  of  the  energy  released  by  the  detona¬ 
tion  of  a  nuclear  weapon  in  the  atmosphere  will  appear  in 
two  primary  forms.  A  blast  wave  will  propagate  outward 
from  the  burst  point  at  supersonic  speeds  and  will  induce  a 
motion  to  the  air  as  well  as  change  the  air  pressure  and 
temperature.  The  hot  fireball  will  radiate  thermal  energy 
which  will  raise  the  temperature  of  the  ground  and  other 
surrounding  surfaces.  The  soil  may  respond  mechanically  to 
this  thermal  disturbance  and  portions  of  it  may  be  thrown 
upward  into  the  surrounding  quiescent  air.  Furthermore,  the 
air  blast -induced  winds  may  scour  the  soil  surface  and 
transport  heated  particulate  matter  upward,  thereby  forming 
a  dust  cloud.  The  rising  motion  of  the  fireball  may  also 
generate  surface  winds,  sometimes  called  "afterwinds",  of 
considerable  magnitude.  These  winds  will  also  contribute  to 
the  generation  or  continued  growth  of  the  dust  cloud.  The 
soil  particulate  matter  sucked  up  by  the  fireball  motion 
will  rise  to  great  heights  and  then  fall  out,  forming  another 
type  of  dust  cloud--one  which  exists  over  great  distances. 
Finally,  for  sufficiently  low  heights  of  bursts,  crater 
ejecta  consisting  of  both  large  and  small  particles  will  be 
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deposited  over  a  considerable  surface  region  surrounding 
the  burst  point. 

The  coarse  part'  ulate  matter  will  settle  out  of 
the  air  rapidly  due  to  gravitational  effects;  however  the 
very  fine  soil  particles  will  remain  airborne  fcr  a  consid¬ 
erable  period  of  time  compared  with  the  duration  of  the 
blast  winds.  These  fine  particles  will  be  transported  hori¬ 
zontally  large  distances  due  to  the  ambient  surface  winds. 

Preliminary  studies'* *  ^  have  indicated  that  the 
major  source  of  dust  at  the  potential  sites  of  interest  will 
be  that  due  to  the  blast -induced  winds.  Thus,  a  theoretical 
model  of  the  dust  cloud  generated  by  the  blast -induced  winds 
has  been  developed  and  is  discussed  in  this  chapter.  The  in¬ 
fluence  of  the  thermal  effects  has  been  neglected  in  this 
dust  model.  However,  the  influence  of  the  heated  dust  and 
other  thermal  effects  is  treated  in  a  more  complex  model 
used  to  predict  both  the  air  temperature  and  dust  density  at 
the  point  of  interest.  The  simplified  or  cold  dust  model  is 
useful  for  generating  data  for  intermediate,  as  well  as  late, 
times  when  the  thermal  effects  become  insignificant.  This 
model  also  enables  one  to  evaluate  the  relative  influence  of 
certain  assumptions,  variables,  and  coefficients  associated 
with  the  dust  cloud. 

The  cold  dust  problem  as  well  as  the  overall  dust 
and  thermal  problem  is  extremely  complex.  It  is  clear  that 
if  engineering  type  estimates  are  to  be  made  of  the  local 
environment,  a  number  of  simplifying  assumptions  must  be  made. 
These  assumptions,  hopefully,  must  be  such  that  the  essential 
of  predominant  features  of  the  process  are  retained  as  well 
as  introducing  a  reasonable,  but  not  excessive,  degree  of 
conservatism  to  the  predictions  or  estimates. 
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Ine  air  blast  environment  in  the  absence  of  the  dust 
and  thermal  effects  is  in  itself  a  relatively  complex  problem, 
involving  the  nonlinear  behavior  of  air  in  a  transient-multi¬ 
dimensional  space.  A  description  of  the  air  pressure,  air 
»  temperature,  and  particle  velocity,  in  some  relatively  simple 

burst  conditions  and  associated  geometric  configurations,  has 
been  established,  using  theoretical  considerations  as  well  as 
experimental  observations  made  during  field  tests  cf  nominal 
size  nuclear  weapons.  The  inclusion  of  the  dust  and  thermal 
aspects  into  the  nonlinear  gas  dynamic  problem  would  render 
the  problem  unsolvable  within  the  allowable  effort.  However, 
reasonable  estimates  ^.f  the  dust  environment  can  be  made, 
provided  it  can  be  assumed  that  the  gas  dynamic  fields  are 
essentially  the  same  as  those  existing  in  the  absence  of  the 
dust.  Such  an  uncoupling  assumption  is  made  herein  and  is 
valid  whenever  the  kinetic  energy  of  the  dust  is  small  compared 
with  the  total  energy  locally.  Thus,  the  particle  velocity 
field  and  pressure  field  are  assumed  locally  and  are  those 
established  from  the  usual  weapons  effects  data.  The  vertical 
transport  of  the  dust  is  assumed  to  be  due  to  the  turbulent 
action  of  the  air;  therefore,  it  is  assumed  herein  that,  at 
least  locally,  the  motion  of  the  air  is  relatively  uniform 
both  in  its  gross  motion  and  in  its  turbulent  intensity. 

This  assumption  implies  that  we  are  dealing  with  flows  near 
relatively  flat  surfaces.  The  model  will  tend  to  break  down 
in  the  neighborhood  of  large  surface  perturbations  which  would 
give  rise  to  strong  local  effects.  Since  we  are  concerned 
primarily  with  local  effects,  we  can  neglect  the  small  hori¬ 
zontal  and  vertical  gradients  which  exist  in  the  gas  dynamic 
fields . 

Section  5.2  presents  the  details  of  the  IITRI  dust 
model  and  a  similarity  solution  for  a  uniform  soil  condition. 
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Another  section  (Section  5.3)  presents  an  evaluation  of  the 
coefficients  used  in  the  dust  model,  together  with  a  comparison 
of  some  limited  field  data.  This  is  followed  by  Section  5.4 
which  presents  a  series  of  typical  results  illustrating  the 
influence  of  the  many  variables  which  contribute  to  the 
characteristics  of  the  dust  cloud.  A  list  of  symbols  used 
in  the  dust  model  is  given  in  Section  5.5  and  Section  5.6  is 
a  brief  estimate  cf  the  dust  cloud  generated  by  crater  ejecta. 

5.2  TITRI  Dust  Cloud  Model 

5.2.1  General 

The  horizontal  flow  of  air  over  the  soil  surface 
will  result  in  aerodynamic  forces  acting  upon  portions  of  the 
soil  mass  or  upon  individual  soil  particles.  If  the  forces 
are  of  sufficient  magnitude,  the  soil  particle  or  soil  mass 
will  move  from  its  original  position  and  bounce,  roll,  or 
otherwise  be  transported  in  the  direction  of  the  flow.  During 
this  erosion  process  the  soil  particles  may  obtain  a  vertical 
component  of  motion  and  thus  become,  at  least  momentarily, 
airborne.  The  cloud  will,  in  addition  to  its  grossly  horizon¬ 
tal  motion,  be  characterized  by  some  degree  of  local  turbulent 
motion.  The  soil  particles,  once  airborne,  can  then  be 
transported  horizontally  by  the  gross  motion  of  the  air  and 
vertically  by  the  inherent  turbulent  motion.  Thus,  one  can 
expect  that  the  blast-induced  winds  will  erode  some  soil 
surfaces  and  transport  the  soil  particles  into  the  air, 
forming  a  dust  cloud.  As  the  intensity  of  the  wind  and  the 
turbulence  decay,  gravitational  forces  will  then  permit  the 
particulate  matter  to  settle  out  of  the  cloud. 

The  very  fine  soil  particles  will  remain  airborne 
for  periods  of  time  which  are  much  greater  than  the  durations 
of  the  blast-induced  winds.  These  particles  will  be  trans¬ 
ported  great  distances  by  the  ambient  surface  winds  which  will 
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generally  exist,  hence  the  dust  environment  at  a  particular 
location  may  be  influenced  by  the  dust  clouds  generated  at 
some  distant  locations. 

"Soil"  is  a  term  used  to  describe  a  large  variety 
of  complex  and  varied  substances.  The  soil  can  be  charac¬ 
terized  by  its  physical  appearance,  physical  properties, 
e hemic 2!  composition,  geological  origin,  etc.  and  by  the 
mass  distribution  of  these  quantities.  Also,  certain  time- 
dependent  characteristics  such  as  moisture  content  and  com¬ 
paction  may  be  exibited.  Furthermore,  the  soil  characteris¬ 
tics  may  change  drastically  with  depth  near  the  soil  surface. 
Since  the  soil  is  such  a  complex  entity,  we  must  describe  it 
in  a  few  simple  terms  and  adjust  the  dust  cloud  predictions 
for  some  of  the  other  soil  characteristics  in  a  gross  manner. 
We  will  therefore  describe  the  soil  as  a  particulate  matter 
with  a  certain  discrete  particle-size  distribution,  particle 
density,  and  bulk  or  in  situ  soil  density. 

The  following  two  basic  assumptions  are  made  in 
the  IITRI  dust  cloud  model.  First,  it  is  assumed  that  the 
horizontal  motion  of  the  dust,  once  ejected  into  the  air,  is 
identical  to  the  gross  horizontal  motion  of  the  air.  This 
assumption  is  valid  for  a  wide  range  of  particle  sizes  pro¬ 
vided  that  the  time  and  space  resolutions  (required  to  ob¬ 
tain  equilibrium  with  the  air  motion)  are  small.  A  particle 
at  rest,  but  in  the  airstream,  will  be  subjected  initially 
to  a  large  drag  force.  The  particle  will  accelerate  rapidly 
until  it  reaches  (asymptotically)  the  air  velocity.  The 
response  time  will  be  of  the  order  of  milliseconds  and  will 
result  in  a  lag  of  approximately  several  feet.  The  second 
basic  assumption  is  that  the  motion  of  any  particle  or  group 
of  particles  is  independent  of  the  existence  or  motion  of 
any  other  particle  or  groups  of  particles.  This  assumption 
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will  be  valid  provided  the  number  of  particles  per  unit 
volume  is  not  excessive  as  is  anticipated.  This  assumption 
permits  one  to  treat  each  group  of  particles  separately  and 
to  obtain  the  total  effect  by  summing  all  of  the  groups. 

5.2.2  The  Erosion  Process 

The  first  phase  of  the  dust  cloud  problem  is  the 
dust  pickup  of  erosion  phase.  It  is  postulated  that  those 
soil  particles  which  are  at  the  surface  will  be  subjected 
to  an  aerodynamic  force  acting  primarily  in  the  horizontal 
direction.  These  particles  will  move  horizontally  and  be 
lifted  into  the  airstream.  It  is  clear  that  successive 
layers  of  particles  cannot  move  into  the  airstream  until  the 
preceding  particles  have  moved  some  finite  distance  away  from 
their  original  resting  place. 

The  force  which  acts  on  the  particles  will  be  equal 

to  or  be  some  fraction  of  the  aerodynamic  drag  caused  by  the 

o 

high  wind  velocity.  This  force  will  be  proportional  to  V  x  D 

where  D  is  the  diameter  of  the  particle.  The  acceleration 

o 

will  be  proportional  to  V  / D  due  to  the  mass  of  the  particle, 
and  the  time  for  the  particle  to  move  a  given  distance  (in 
terms  of  its  size  D)  will  be  proportional  to  D/V.  The  number 
of  available  particles  per  layer  will  be  proportional  to 
1/D".  Therefore,  this  is  the  number  of  particles  that  can  be 
stripped  away  or  eroded  in  the  above  time  interval;  hence 


where 


dn 

3t 


(5.1) 


v  =  number  of  particles  available  per  unit  volume, 

§ 

V  =  air  velocity, 

n  =  number  of  particles  in  a  given  size  class,  and 
t  =  time. 


This  relationship  can  be  rewritten  as 


3F  -  apgfp  |V| 


(5.2) 


and 


5vt\  r  5vt" 

|v|  j  [  |V|j* 


(5.3) 


where 


f  =  erosion  cutoff  criteria, 
a  =  erosion  constant, 
w  =  weight  of  particle  per  unit  area, 
os  =  specific  weight  of  soil, 
t  -  t ime , 

p  =  weight  fraction  of  particles, 

V  =  air  velocity, 

Vt  =  terminal  velocity  of  particle,  and 
)-( ,  ;  =  unit  step  function. 

_u  J 

The  factor  f  (used  in  Equation  5.2)  represents  a 

f 

cutoff  mechanism  which  states  that  the  air  velocity  must  be 

some  multiple  of  the  terminal  velocity  (5,  in  this  case) 

before  a  particle  can  leave  the  surface.  Some  experimental 
5  2 

evidence  '  exists  to  support  such  a  concept,  at  least  for 
loose-sand  surfaces.  The  cutoff  function  is  arbitrarily  made 
smooth  by  multiplying  the  unit  function  by  the  factor 
(1  -  5Vt/|V|). 

The  factor  a  is  the  erosion  constant  and  has  been 
evaluated  for  some  soils  during  the  current  program.  These 
experiments  are  described  in  Appendix  A.  A  value  of  6  x  10  ^ 
has  been  selected  as  being  applicable  to  the  current  problem. 


5.2.3 


The  Transport  Process 


The  transport  of  the  dust  particles  is  postulated 
to  be  due  to  gross  horizontal  motion  of  the  air,  as  well  as 


125 


being  due  to  local  turbulent  motion.  It  is  possible  now  to 
consider  the  motion  of  a  small  group  of  particles  (of  the 
same  size)  which  are  ejected  into  the  air  at  a  range  of 
R  =  and  a  time  of  t  =  t^  (see  Figure  5.1).  The  weight  of 
this  group  of  particles  is  dw.  According  to  the  first  basic 
assumption,  the  horizontal  motion  of  this  group  of  particles 
is  identical  to  the  horizontal  motion  of  the  air  (identified 
in  Figure  5.1  as  air  column  A).  The  air  velocity  is  assumed 
uniform  in  the  vertical  direction,  and  no  boundary  layer 
effects  are  included.  This  assumption  implies  that  the 
heights  of  interest  are  generally  much  greater  than  any 
boundary  layer  thickness.  The  individual  particles  of  this 
group  or  dust  packet  will  also  move  in  the  vertical  direc¬ 
tion  due  to  the  action  of  air  turbulence  and  gravitational 
forces.  However,  the  vertical  motion  will  be  bounded  by  some 
cloud  height  H.  This  is  illustrated  in  Figure  5.1  at  t  s  t2 
and  t  =  t^.  The  rate  at  which  the  cloud  height  grows  (in 
the  absence  of  gravity)  will  be  approximately  proportional  to 
the  absolute  magnitude  of  the  gross  (horizontal)  air  velocity, 
'"’he  cloud,  and  therefore  its  extremity  H,  will  also  settle 
uva  to  the  action  of  gravity.  This  decreasing  rate  will  be 
equal  to  the  terminal  velocity  Vt  of  the  particle.  Thus, 
we  conclude  that  the  time  history  of  the  cloud  height  as¬ 
sociated  with  the  particular  dust  packet  is  given  by 

HI  =  K|V|  -  Vt  .  (5.4) 

where  V  is  the  instantaneous  horizontal  wind  velocity  of 
air  column  A  and  K  is  a  factor  yet  to  be  evaluated.  The 
value  of  K  will  be  some  small  fraction  of  unity  since  the 
turbulent  velocity  component  will  be  some  small  fraction  of 
the  mean  air  velocity  V.  One  can  also  postulate  that  the 
factor  5  used  in  Equation  5.3  is  related  to  the  factor  K,  as 
K  =  1/2.  This  implies  that  the  erosion  cutoff  depends  upon 
the  turbulent  intensity  of  the  air. 
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Air  Column 


The  distribution  of  the  particles  within  the  cloud 
is  needed  to  determine  the  number  of  particles  per  unit  volume 
(i.e.,  the  density)  at  any  height  y.  It  is  assumed  that  the 
normalized  distribution  (or  probability  function)  is  similar 
at  all  times  and  for  all  particle  sizes.  At  the  present  time 
there  is  little  evidence  to  support  the  selection  of  any 
particular  distribution  function;  however,  one  can  make  an 
intuitive  selection  and  then  evaluate  the  sensitivity  of  the 
results  to  this  selection.  The  simplest  distribution  func¬ 
tion  that  one  might  choose  is  a  uniform  (rectangular)  density 
variation.  Such  a  distribution  is  unlikely,  and  it  is  proba¬ 
ble  that  most  of  the  dust  particles  will  concentrate  at  the 
lower  end  of  the  cloud.  It  is  expected  that  the  distribution 
will  not  be  irregular,  nor  that  the  density  will  be  vanishingly 
low  near  the  upper  portions  of  the  cloud.  If  the  latter  were 
the  case  then  one  could  re-evaluate  the  constant  K  and  select 
a  new  e: .ective  height. 

The  next  simplest  distribution  selected  was  the 
triangular  distribution  in  which  the  density  decreases  line¬ 
arly  with  increasing  height,  assuming  a  value  of  zero  at  the 
top  of  the  cloud.  Thus,  the  density  dp  is  given  by 

dp  =  dp  (0)  1  -  d  0  <  y  <;  H,  (5.5a) 

and 

do  =  0  H  <  y,  (5.5b) 

where  do (0)  is  the  value  of  do  at  y=0.  Since  the  total 
weight  dw  of  the  dust  packet  is  known,  do (0)  can  then  be 
evaluated.  This  yields 

do  =  —jj-*  1-^  0<:y<H.  (5.6) 
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The  distribution  can  be  generalized  somewhat  by  assuming  that 


do  =  do  (0) 


1 


Or  y  -  H. 


(5.7) 


However,  this  introduces  another  unknown  parameter  c,  which 
must  be  evaluated.  Until  some  experimental  evidence  is 
developed,  the  triangular  form  (c=0)  will  be  used.  In  a  later 
part  of  this  section  a  series  of  idealized  problems  for  each 
of  the  above  distribution  functions  is  evaluated  showing  that 
the  results  are  not  significantly  influenced  by  the  selection 
of  the  distribution  function.  This  is  particularly  true 
between  the  triangular  form  (c=0)  and  the  concave  triangle 
(c=l),  the  more  realistic  distribution  functions. 

The  terminal  velocity  plays  an  important  role  in 
these  dust  studies.  For  particle  sizes  greater  than  0.1  mm 
(100  microns)  in  diameter,  Stokes'  equation  is  inadequate. 
Stokes'  equation  is  valid  for  Reynolds  number  N  of  unity  or 
less.  The  drag  coefficient  Cj  for  spheres  is  presented  in 
Figure  5.2.  The  drag  coefficient  decreases  rapidly  with 
increasing  Reynolds  number  and  reaches  a  constant  value  in 

the  turbulent  flow  region.  We  have  approximated  the  experi- 

5  3  . 

mental  results  *  by  an  equation  of  the  form 


r  _  24  ,  1 

Ld  "  TT  '  7  * 

e 


(5.8) 


This  simple  form  appears  adequate  over  the  large  Reynolds 
number  range  of  concern  (i.e.,  corresponding  to  particles 
up  to  the  size  of  medium  gravel).  The  terminal  velocity  Vt 
is  given  by 
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(5.9) 


1 


where 

D  *  diameter  of  particle, 
li  *  viscosity  of  air, 
oa  =  density  of  air,  and 
o  ,  =  density  of  particle. 


5.2.4  Blast-Induced  Wind 


The  blast -induced  winds  are  a  function  of  the  many 
weapons  effects  variables.  The  local  wind  variation  or  velo¬ 
city  field  can  be  expressed  in  terms  of  a  number  of  coeffi¬ 
cients,  which  in  turn  must  be  evaluated  from  the  specific 
nuclear  attack  situation.  For  a  simple  case,  the  velocity 
variation  Vr  of  the  air  at  a  fixed  range  is  generally  assumed 
to  be  of  the  form 


- » 
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I-  kt 

t 

L  °J 

O 

where 

VQ  *  peak  air  velocity, 
t  *  positive  phase  duration,  and 
k  *  exponential  coefficient. 


(5.10) 


Generally,  this  equation  is  valid  well  into  the  negative 
phase.  Furthermore,  it  can  be  demonstrated  that  the  maximum 
horizontal  displacement  of  the  air  at  the  low  overpressure 
range  will  be  of  the  order  of  one -sixth  of  the  slant  range; 
and  that  over  such  an  increment  of  range,  the  velocity  field 
will  not  vary  appreciably  in  the  horizontal  direction.  There 
fore,  one  can  construct  an  expression  for  the  velocity  in 
terms  of  range  and  time  from  Equation  5.10  and  from  an  assump 
tion  that  the  field  is  composed  of  simple  waves.  For  the 
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current  calculations,  however,  we  need  to  know  the  velocity- 
time  variation  along  an  air  particle  path.  Its  derivation 
from  the  above  indicated  velocity  field  would  be  somewhat 
cumbersome.  To  avoid  this  complication,  it  is  assumed  that 
along  an  air  particle  path  the  velocity  is  of  the  form 


(5.11) 


where 

to 


is  a  modified  positive  phase  duration  and  is  equal 


Ro 

i r 


(5.12) 


where 

R  =  maximum  displacement  of  air,  and 
U  =  air  shock  velocity. 

The  positive  phase  duration  has  been  increased,  to  account 
for  the  additional  time  that  the  particle  is  in  motion.  Fig¬ 
ure  5.3  illustrates  the  air  particle  path  and  shows  that  the 
air  particle  originally  at  is  displaced  a  maximum  distance 
Rq  and  reaches  this  position  after  moving  for  a  time  equal 

to  the  sum  t,  +  t  .  The  time  t,  is  the  transit  time  of  the 
1  o  1 

shock  front  over  the  range  increment  Rq.  The  range  R*RQ  is 
the  range  of  interest.  The  problem  is  to  determine  the  dust 
history  at  this  particular  position.  The  dust  that  exists 
at  the  location  of  interest  is  that  dust  which  is  picked  up 
along  the  particular  air  particle  path  which  intersects  RQ. 

A  simple  integration  of  all  of  the  elemental  contributions 
along  these  paths  yields 


P 


2 

H 


aPgfp|v| dt 


0  s  Y  *  H,  (5.13) 
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and 

p  =  0  H  *  Y,  (5.14) 

where 

Y  =  height  of  interest,  and 
ta  =  time  of  arrival  of  air  shock  at  R^. 

This  integration  is  carried  out  for  each  particle  size  group. 
It  should  be  noted  (Figure  5.3)  that  some  air  particle  paths 
(between  Path  A  and  Path  C)  cross  the  plant  location  twice, 
whereas  others  only  cross  the  plant  location  once.  The  late¬ 
time  dust  conditions  at  the  plant  are  those  accumulated  along 
Path  C.  Finally,  one  can  place  limits  on  the  range  incre¬ 
ments  from  which  the  dust  source  may  exist  and  still  influence 
the  dust  cloud  at  R=RQ.  The  range  increment  RQ  is  given  by 

V  t* 

R0  =  -2j2  k  -  1  +  e-kJ  ,  (5.15) 

and  the  range  R^  is  given  by 

V  t* 

R4  =  2R0  -  -2^2  (k  -  1)  .  (5.16) 

The  range  increment  R^  for  a  wide  range  of  weapon  corditions 
is  of  the  order  of  3000  ft  or  less. 

Once  the  air  velocity  has  been  established  we  can 
integrate  Equation  5.4  and  determine  the  cloud  height  of  any 
particular  packet  of  dust;  namely, 

H(t1,t,Vt)  -  H*(t)  -  H*(tt)  -  Vt(t  -  tx)  (5.17) 


134 
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where 


H*(t)  = 


KV  t* 
o  o 


l 


(k  -  1) 


(5.18) 


when  t  <  t  ;  and 
o 


KV  t* 

_  ovd  _  o 


H*(t)  =  2KRo  - 


k2 


-  1)  + 


kt 

L  O 


(5.19) 

when  t*  £  t.  A  typical  example  is  presented  in  Figure  5.4 
which  illustrates  the  effect  of  particle  size  on  the  height  of 
the  cloud.  It  should  be  noted  that  the  cloud  height  decreases 
and  vanishes  quite  rapidly  for  the  heavier  particles. 


5.2.5  Analytical  Solutions 

A  series  of  analytical  solutions  can  be  readily 
obtained  for  the  case  of  uniform  ground  conditions  where  the 
soil  properties  do  not  vary  in  the  horizontal  direction. 

These  solutions  are  useful  in  examining  the  influence  of  some 
of  the  many  variables,  and  in  some  instances  will  be  useful 
for  site  evaluations.  In  particular,  the  sensitivity  of  the 
dust  density  to  the  form  of  the  distribution  can  be  examined. 

The  assumption  of  uniform  ground  conditions  elimin¬ 
ates  the  need  of  knowing  the  position  of  the  air  column  at 
any  stage  of  the  flow,  and  indicates  that  all  of  the  air 
columns  behave  identically.  We  can  obtain  solutions  for  each 
particle  class  separately  (let  p=l)  and  combine  the  individual 
for  any  particle  size  distribution  of  interest. 

A  general  solution  can  be  obtained  by  recalling 

aPsp 

dw  =  (K|V|  -Vt)  dt,  (5.20) 


solutions 

that 
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and  using 


(5.21) 


to  obtain 


aP  sp 

dw  —  dH  . 


(5.22) 


The  elimination  of  the  specific  function  defining  the  air 
velocity  is  the  key  to  this  simplification  and  can  be  accom¬ 
plished  by  equating  the  factor  5  (Equation  5.3)  to  1/K  as 
discussed  in  the  previous  subsection.  Thus,  Equation  5.22 
can  be  integrated  over  all  dust  packets  to  obtain  the  density 
of  a  given  size  class. 


(5.23) 


where  the  limits  of  integration  for  H  are  from  the  top  of 
the  cloud  H*  to  the  elevation  Y  of  interest.  These  integra¬ 
tion  limits  are  illustrated  in  Figure  5.5.  The  integration 
can  be  carried  out  over  the  variable  t^  or  over  the  variable 
H.  The  integration  of  Equation  5.23  yields  a  similarity  solu 
tion  in  the  variable 


that  is , 

P  (C) 


2ao_p 

~-(C  -  1  -  In  ;) 


(5.24) 

(5.25) 


H*  is  a  function  of  the  time  measured  along  a  particle  path. 
To  convert  to  an  observation  time  t 1  at  a  fixed  position,  one 
must  use  the  relation 
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t 

t'  =  t  -  -J  f  V(t)  dt  .  (5.26) 

J0 

H*  is  obtained  by  integrating  Equation  5.2  as 

t 

H*  =  f  (Kj V J  -V.  )  dt  .  (5.27) 

t) 


For  the  assumed  velocity  field  discussed  previously,  these 
become 


(5.28) 


where  HV (t )  is  defined  by  Equations  5.16  and  5.17.  The  above 
similarity  solution  is  based  on  the  triangular  distribution 
function.  Corresponding  solutions  have  been  obtained  for 
both  the  rectangular  and  concave  triangle  distribution  func¬ 
tions,  respectively,  as 


a;:,sP 

p(Y,t)  =  — £ — In  r 


(5.29) 


and 

o(Y,t) 


r 


=  B 


-1 


-  e”“  -  In  C  + 


a  n 


-  1X--1) 

_  _  i 


n-1 


n.n. 


n=i 


(5.30) 


where 


The  dust  density  variation  with  height 


(5.31) 

(actually  the  similar¬ 


ity  variable  C)  is  presented  in  Figure  5.6  for  each  of  the 
assumed  distribution  functions.  This  variation  clearly  indi¬ 
cates  that  the  density  does  not  differ  appreciably  between 
the  triangular  and  the  concave  triangle  (c=l)  distribution 
functions  except  at  extremely  low  heights.  These  low  heights 
are  of  no  interest  to  the  current  problem  and  represent  the 
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dust  boundary  layer  immediately  adjacent  to  the  ground  surfact. 
The  influence  of  the  distribution  function  on  the  time  details 
of  the  density  is  presented  in  Figure  5.7  for  a  representa¬ 
tive  set  of  parameter  values.  The  results  indicate  that  for 
multisecond  time  intervals,  the  triangular  distribution  is 
more  than  adequate . 

5.2.6  Late  Time  Dust  Cloud  Model 

The  blast  winds  which  accompany  the  detonation  of 
a  nuclear  weapon  in  the  atmosphere  will  generate  wind  induced 
dust  clouds  which  cover  a  very  large  area  and  which  will 
remain  airborne  for  a  substantial  period  of  time.  The  proper¬ 
ties  of  the  dust  cloud  will  vary  throughout  the  cloud  forma¬ 
tion  due  to  the  local  variations  of  the  surface  soil  condi¬ 
tions  and  the  variation  in  the  magnitude  and  duration  of  the 
blast  winds.  The  entire  cloud  formation  will  also  be  trans¬ 
ported  over  the  terrain  due  to  the  ambient  surface  winds. 

The  larger  soil  particles  which  are  initially  raised  by  the 
strong  blast  winds  will  settle  out  of  the  cloud  at  the  latter 
stages  when  the  blast  winds  subside.  However,  the  very  fine 
soil  narticles  will  remain  airborne  and  will  be  further 
dispersed  due  to  the  ambient  winds.  All  of  the  airborne 
soil  particles  will  be  transported  in  a  gross  manner  hy  the 
surface  wind  and  will  be  accompanied  by  simultaneous  turbulent 
diffusion  in  all  directions.  Other  mechanisms  will  also 
contribute  to  the  further  local  transport  and/or  mixing  of 
the  dust  particles.  These  other  mechanisms  would  include 
the  influence  of  local  thermal  effects,  the  blast  after 
winds  when  they  occur,  and  other  local  meteorological  phenome¬ 
non  such  as  overturning. 

The  late-time  dust  effects  are  important  since  much 
of  the  mechanical  and  thermal  air  breathing  equipment  used 
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in  hardened  facilities  must  operate  during  the  post-attack 
period.  The  fine  details  of  the  dust  cloud  are  probably  not 
important  for  long  time  performance  of  this  type  of  equip¬ 
ment.  Rather,  average  values  of  the  dust  density  variations 
should  be  sufficient.  Thus,  this  complex  problem  can  be 
simplified  considerably.  A  late-time  dust  model  has  been 
written  and  used  to  examine  a  variety  of  typical  attacks  and 
meterological  conditions.  This  late-time  dust  cloud  model 
utilizes  an  average  or  uniform  ground  condition  for  the  area 
covered  by  the  entire  cloud  formation.  The  ambient  surface 
winds  are  assumed  to  be  uniform  in  both  magnitude  and  direc¬ 
tion  for  the  period  of  interest.  Horizontal  dispersion,  as 
well  as  other  transport  and  mixing  mechanisms,  have  been  ex¬ 
cluded  in  this  model.  The  neglect  of  the  horizontal  disper¬ 
sion  is  reasonable  on  the  basis  that  for  uniform  ground  condi¬ 
tions,  the  horizontal  gradients  of  the  dust  concentrations 
throughout  the  entire  cloud  formation  will  not  be  large. 

The  vertical  gradients  are  at  least  an  order  of  magnitude 
greater  than  the  horizontal  gradients.  The  neglect  of  other 
transport  and  mixing  mechanisms  will  result  in  a  conservative 
estimate  of  the  dust  environment  that  is  in  an  upper  bound 
value.  For  low  altitude  burst  conditions,  crater  ejecta, 
fireball  dust,  and  dust  raised  by  afterwinds  will  influence 
the  dust  environment,  the  latter  two  being  late-time  effects. 
The  soil  trapped  in  the  fireball  will  be  carried  off  and  most 
likely  be  transported  a  great  distance  from  the  burst  area 
before  the  soil  particles  settle  out.  All  of  these  effects 
have  been  excluded  from  the  present  model. 

The  early  phase  of  the  growth  and  the  settling  of 
the  dust  cloud  at  any  arbitrary  position,  such  as  at  the  site, 
has  been  determined.  Typically  the  shock  front  emanating 
from  the  detonation  will  arrive  at  the  site  at  times  of  the 
order  of  tens  of  seconds.  The  positive  phase  duration  of  the 
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blast  wind  will  be  of  the  order  of  many  seconds  such  that  the 
primary  growth  or  establishment  of  the  entire  dust  cloud 
formation  will  be  completed  in  a  very  short  time  period. 

Thus,  dust  densities  at  a  site  typically  reach  their  maximum 
values  prior  to  50-100  seconds  after  detonation.  Thus,  for 
late-time  information,  one  can  assume  that  the  entire  dust 
cloud  formation  is  established  instantaneously  and  at  the 
time  of  detonation.  Figure  5.8  illustrates  the  dust  cloud 
formation  and  identifies  the  position  of  the  site  relative 
tc  the  burst  point  or  "ground  zero"  (the  origin).  The  dust 
cloud  formation  will  be  symmetrical  about  the  vertical  axis 
through  the  burst  point. 

The  ambient  surface  winds  are  assumed  to  be  uniform 
in  both  magnitude  and  direction  such  that  the  cloud  formation 
will  drift  past  the  site  and  create  an  effective  path  of  the 
site  relative  to  the  cloud  formation.  This  straight  line 
path  is  illustrated  in  Figure  5.8(a)  and  is  defined  by  the 
angle  0  measured  between  the  negative  wind  vector  and  a  line 
generated  by  the  site  location  and  ground  zero.  The  distance 
Rt  measured  along  this  path  is  given  by 

Rt  =  W  •  t  (5.32) 

where 

W  =  wind  velocity  and 

t  =  time. 

The  range  R  or  distance  from  the  origin  of  the  cloud  forma¬ 
tion  of  a  point  along  this  path  is  given  by 

R  ®*VRt  +  Rs  ’  2RsRt  cos  0  (5.33) 

where 

R  =  initial  range  of  the  site  (t  =  0). 
s 
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(a)  Plan  View 


Figure  5.8  Dust  Cloud  Formation 
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The  characteristics  of  the  dust  cloud  are  dependent 
upon  both  the  location  R  within  the  dust  cloud  and  the  total 
elapsed  time  t.  The  location  defines  the  overpressure  and 
hence,  the  local  blast  wind  parameters  which  generate  the 
dust  cloud;  the  total  elapsed  time  defines  the  degree  of 
cloud  growth  or  settling  which  occurs  at  late  times. 

The  previous  subsections  described  a  dust  model  for 
the  early  phases  of  the  dust  cloud  formation.  The  similarity 
solution  for  a  uniform  ground  condition  will  be  used  in  this 
late-time  dust  model;  however,  some  modifications  will  be 
introduced  in  order  to  better  describe  the  late -time  dust 
environment.  The  late -time  model  is  restricted  to  a  single 
burst  situation  and  to  a  uniform  ground  condition.  Neither  of 
these  restrictions  is  absolute;  however,  their  adoption 
greatly  simplifies  the  model  and  the  subsequent  computational 
effort  involved  in  evaluating  a  particular  problem.  The 
similarity  solution  and  its  subsequent  modifications  apply  to 
each  specific  particle-size  class  of  interest  or  those  avail¬ 
able  in  the  soil.  The  total  dust  density  is  equal  to  the  sum 
of  the  dust  density  contributions  from  each  of  the  individual 
particle  size  classes  used. 

The  magnitude  of  the  ambient  surface  winds  will 
generally  be  sufficiently  low  such  that  additional  surface 
soil  erosion  or  soil  "pickup”  will  not  occur.  It  is  there¬ 
fore  assumed  that  no  soil  particulate  matter  enters  the  cloud 
after  the  actual  blast  wind  phase  of  the  cloud  development. 

The  initial  conditions  or  state  of  the  dust  cloud  is  therefore 
defined  by  the  blast  wind  parameters.  The  total  amount  of 
dust  for  each  size  class  and  the  initial  height  of  the  cloud 
formation,  together  with  the  distribution  of  the  particulate 
matter,  are  sufficient  to  define  the  initial  state  of  the  dust 
environment . 
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The  initial  height  of  the  dust  cloud,  neglecting 
settling,  is  given  by  the  limiting  value  of  H^(t)  as  t  -*■  *. 
Equation  5.19  yields 


lim  H*(t)  =  2KjR0  -  -  ;°  (k  -  1)  =  H*  . 
«d  k 


(5.34) 


The  initial  dust  density  for  each  particle  size  class  is 
given  by 

2a°cP  „  / 


where 


oC)  = 


c  -  y/h;  . 


(S  -  1  -in£) 


(5.35) 


(5.36) 


The  rate  of  change  of  the  cloud  height  as  given  by 
Equation  5.4  is  still  applicable;  however,  the  ambient  surface 
wind  W  is  now  the  driving  force.  We  therefore  obtain 


a?  - ««  -  vt  • 


(5.37) 


If  dH/dt  is  positive,  then  further  growth  of  the  cloud  (i.e., 
for  the  specific  particle  size  class  associated  with  the  value 
of  V  )  will  occur.  Since  no  further  erosion  has  occurred 
since  the  initial  cloud  was  formed,  the  similarity  solution 
does  not  apply  directly.  At  some  time,  t,  the  cloud  will 
reach  a  new  height,  H,  given  by 


where 


H  -  H0  -  Ct 


C  *  Vt  -  KW  . 


(5.38) 


(5.39) 


It  is  assumed  that  the  particulate  matter  is  distributed  over 
the  height  H  in  the  way  that  it  was  distributed  originally 
over  the  height  H*,  but  conserving  the  overall  mass  of 
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the  cloud.  The  dust  density  for  the  growth  phase  (C  <  0) 
is  given  by 

* 

*  H  A 

0  =  0  TT  Cc  '  1  *^n  c]  (5.40) 

whc’*e 

•jjf 

a  =  (2ao  p)/K  and 
C  =  Y/H. 


For  comparison  purposes  it  will  be  convenient  to  recast  this 
equation  in  terms  of  a  dimensionless  time  parame^fc,  T. 
Therefore , 


C, 


c  = 


rr^rr 


(i  + 1  > 


-  1  */n 


where 


and 


Co  -  r*  >  0  s  Co  s  1 

O 


Ct 


H 


,  0  £  T  . 


,  (5.41) 


(5.42) 


(5.43) 


o 


The  density  distribution  for  the  growth  phase  is 
illustrated  in  Figure  5.9(a)  and  indicates  that  the  density 
in  the  lower  portion  of  the  cloud  decreases  with  time  as  the 
particulate  matter  migrates  upward.  The  dust  density  in  the 
upper  portion  of  the  cloud  will,  on  the  other  hand,  increase 
with  time.  The  oust  density  at  some  height  which  is  greater 
than  the  initial  height  of  the  cloud  will  be  zero  until  the 
cloud  heig  »t  grows  sufficiently  to  reach  this  elevation. 
Generally,  ^nlerest  will  be  limited  to  the  lower  portions  of 
the  initial  cloud.  The  dust  density  variation  with  time  for 
the  growth  phase  is  illustrated  in  Figure  5.10  for  several 
values  of  £  .  The  dust  density  varies  quite  slowly  with  time. 

V 

The  cloud  height  has  doubled  for  a  value  of  t  *  1. 
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Figure  5,9  Late  Time  Dust  Cloud  Models 
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Figure  5.10  Late  Time  Dust  Cloud  Density  Models 
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If  dH/dt  is  negative  (C  >  0),  then  the  cloud  height 
will  decrease  and  settling  of  the  particulate  material  will 
occur.  In  this  case,  the  similarity  solution  can  be  used,  viz. 

c=osL'“l-£n']  (5.44) 


where 


(5.45) 


In  terms  of  the  dimensionless  time,  * , 
becomes 

,1 

0  =  0  ~  / 

where 


-o 

T77 


-  1  - 


Equation  (5.44) 


(5.46) 


b  o 


(5.47) 


and 

T  =  %  ,  0  :  t  .  (5.48) 

Ho 

Note  that  Equations  (5.43)  and  (5.^d)  differ  in  sign.  This 
difference  is  due  to  desire  to  make  t  increase  with  increas- 
int  t  while  the  sign  of  C  differs  between  the  growth  and 
settling  phase. 

The  above  similarity  settling  equation  describes 
the  fact  that  the  dust  particles  at  the  upper  portions  of  the 
cloud  fall  faster  than  those  near  the  bottom  of  the  cloud. 
Thus  the  dust  begins  to  concentrate  near  the  bottom  of  the 
cloud.  The  dust  density  decreases  with  time  since  some  par¬ 
ticles  do  settle  out  on  the  ground  surface.  It  would  appear, 
however,  that  another,  perhaps  more  realistic  model  could  be 
adopted,  namely  one  in  which  all  of  the  dust  particles  settle 
at  the  same  rate.  This  type  of  behavior  is  referred  to  as 
linear  settling  and  the  dust  density  is  given  by  the  relation 
ship 
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(5.49) 


or  in  terms  of  -  and  '  ,  as 

r  /  v“] 

0  =  °S  |  o  +  1  "  1  ■  Zn  \o  +  T; ;  •  (5.50) 

1_  — 1 

Figure  5.9(b)  and  (c)  illustrates  the  difference  between  simi¬ 
larity  and  linear  settling.  Figure  5.10  presents  the  dust 
history  for  several  values  of  '0  for  each  of  the  settling 
models.  The  linear  settling  model  results  in  a  much  more 
rapid  decrease  in  dust  density  with  time.  In  both  cases  the 
dust  density  vanishes  when  -  =  1  -  ' c. 

Figure  5.11  presents  that  value  of  the  surface  wind 
velocity  for  which  a  particular  size  particle  will  just  settle 
out.  This  value  of  the  wind  velocity  has  been  defined  as  the 
critical  wind  velocity  Wc .  It  should  be  noted  that  particles 
which  are  greater  than  100  to  200  microns  in  size  will  gener¬ 
ally  settle  out  of  the  cloud.  Thus,  the  late  time  dust  envi¬ 
ronment  will  consist  of  only  the  smaller  particles.  Figure  5.12 
presents  a  comparison  between  the  results  of  linear  and  simi¬ 
larity  settling  for  a  typical  problem.  The  discontinuities 
in  the  curves  are  due  to  the  discrete  particle  size  classes 
used  in  the  calculation  and  their  settling  out  at  specific 
times  (r  =  1  -  *0)  for  these  size  classes. 

The  linear  settling  model  described  above  is  the 
settling  model  which  has  been  selected  and  used  in  all  sub¬ 
sequent  calculations. 

5.3  Evaluation  of  Dust  Model  Coefficients 

The  dust  environment  model  is  based  on  a  series  of 
assumptions  utilizing  several  coefficients  to  define  the 
magnitude  of  some  of  the  effects.  It  is  essential  to  examine 
the  validity  of  all  of  the  details  of  the  model  in  order  to 
establish  the  overall  reliability  of  the  prediction. 
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Figure  5,12  Comparison  of  Settling  Models 


Since  the  state  of  the  art  of  generation  and  trans¬ 
port  of  dust  clouds  is  not  sufficiently  advanced,  a  series 
of  two  experiments  has  been  performed  as  part  of  the  current 
program.  One  experiment  deals  with  the  erosion  of  the  soil 
subjected  to  a  transient  surface  wind  and  the  results  will 
assist  in  establishing  the  magnitude  of  the  erosion  constant 
a  used  in  Equation  5.2.  The  second  experiment  deals  with  the 
transport  of  the  dust  cloud,  once  established.  This  experi¬ 
ment  will  determine  the  relationship  between  the  transport 
constant  K  and  the  turbulent  level  of  the  air  stream,  and  will 
also  yield  some  information  on  the  dust  distribution  within 
the  cloud.  Both  experiments  were  rather  limited  in  scope; 
however,  they  are  sufficient  to  substantially  increase  the 
reliability  of  the  prediction  method.  A  description  of  these 
experiments  is  presented  in  Appendix  D.  In  addition  to  these 
experiments,  dust  data  have  been  obtained  from  both  nuclear 
and  high-explosive  field  tests.  These  data,  however,  are  sparse 
and  perhaps  somewhat  unreliable. 

The  erosion  coefficient  a  was  the  most  uncertain 
coefficient  of  the  current  model,  and  the  predicted  dust 
densities  are  directly  proportional  to  its  magnitude.  The 
probable  range  of  its  value,  however,  was  determined.  The 
model  considers  the  magnitude  of  the  aerodynamic  forces  avail¬ 
able  to  remove  layers  of  particles,  and  then  postulates  the 
average  distance  the  particles  must  move  before  the  subsequent 
layers  begin  to  move.  The  estimate  of  the  time  required  to 
remove  each  layer  and  the  number  of  particles  (or  the  weight 
in  each  layer  is  sufficient  to  define  the  magnitude  of 
erosion  constant.  The  following  result  is  obtained  for 
case  where  the  aerodynamic  force  is  constant;  that  is, 
particle  velocity  is  small  compared  to  the  .ir  velocity. 

a  =  V m  (pf)  e  ’  (5'51) 
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whe  re 


fa  =  density  of  the  air, 

j?c  =  density  of  the  particle,  and 

N  =  the  number  of  particle  diameters  traveled 
in  horizontal  direction. 

The  ratio  of  the  air  to  particle  densities  is  well  known  and 
will  not  vary  significantly.  The  relatively  insensitive 
square -root  dependence  of  a  on  N  indicates  that  a  reasonable 
estimate  of  N  will  be  sufficient.  The  value  of  the  erosion 
constant  is  presented  in  Figure  5.13  as  a  function  of  N. 

The  results  of  Equation  5.51  are  indicated  as  Curve  1  in 
this  figure.  Curve  2  is  the  result  of  a  more  refined  analysis 
in  which  the  relative  velocity  between  the  dust  particle  and 
the  air  is  taken  into  consideration  when  reducing  the  eroding 
force . 

A  dust  boundary  layer  will  exist  in  which  the  actu¬ 
al  erosion  process  takes  place.  This  boundary  layer  will  be 
rather  dense  and  the  air  velocity  therein,  which  represents 
the  driving  force,  will  be  substantially  reduced  from  the 
free-field  velocity.  Therefore,  the  erosion  factor  must  be 
reduced  by  the  corresponding  reduction  factor  in  velocity.  The 
dashed  line  in  Figure  5.13  represents  such  a  reduction  corre¬ 
sponding  to  a  factor  of  approximately  0.1. 

The  results  of  the  erosion  experiments  indicated 
that  the  value  of  the  erosion  coefficient  was  considerably 
higher  (approximately  10  )  for  very  loose  soils.  However, 

whenever  a  small  amount  of  compaction  was  applied,  the  value 
of  the  erosion  coefficient  was  reduced  considerably.  The 
results  of  the  erosion  experiments  are  summarized  in  Fig¬ 
ure  5.14.  A  value  of  6  x  10"^  has  been  selected  as  a  reason¬ 
able  value  for  the  erosion  coefficient  for  the  current  prob¬ 
lem,  since  the  surface  soils  will  be  somwhat  compacted. 

The  results  of  the  transport  experiments  indicate 
that  the  value  of  the  transport  coefficients  is  not  strongly 
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dependent  upon  either  the  turbulence  level  or  the  air  velocity. 
These  experimental  results  are  presented  in  Figure  5.15.  A 
value  of  0.2  appears  to  be  a  reasonable  value  for  the  range 
of  variables  and  soil  types  investigated. 

Another  estimate  of  the  value  of  the  transport  coef¬ 
ficient  can  be  obtained  from  a  field  experiment.  The  height 
of  a  dust  cloud  is  one  of  the  few  characteristics  which  can 
be  observed  easily  and  with  reasonable  accuracy.  Observation 
made  during  Shot  7  of  Operation  Greasystake^ is  presented 
in  Figure  5.16,  together  with  a  prediction  based  on  Equation 
5.4.  The  velocity  history  was  estimated  based  on  the  yield, 
the  height  of  burst,  and  the  overpressure  level.  A  value  for 
K  of  approximately  1/4  yields  a  good  comparison  with  the  ex¬ 
perimental  observations.  It  is  significant  that  the  compari¬ 
son  is  good  both  initially  where  the  peak  velocity  is  well 
known  and  at  the  latter  stages  of  the  cloud  growth  (approxi¬ 
mately  one -half  of  the  positive  phase  duration  of  the  veloc¬ 
ity).  This  particular  observation  was  for  a  high-explosive 
test  and  may  not  be  fully  representative  of  a  nuclear 
detonation. 

Several  dust  density  measurements  have  been  made 
of  dust  clouds  raised  by  the  blast  winds  induced  by  low-yield 
nuclear  detonations.  The  reliability  of  these  observations 
is  poor;  however,  these  data  represent  the  only  opportunity 
with  which  to  make  a  comparison  at  this  time  (these  compari¬ 
sons  are  presented  in  Figure  5.17).  These  predictions  were 
based  on  a  uniform  soil  assumption,  and  a  single  particle 
size  (100  microns)  was  used.  It  should  be  noted  that  the 
observations  were  made  during  the  first  half  of  the  positive 
phase  duration  and  no  late -time  effects  were  observed. 

The  experimental  observations  exhibit  rather  large 
local  fluctuation  which  may  be  due  to  either  the  variability 
of  the  soil  in  front  of  the  dust  probe'  or  to  the  interaction 
of  relatively  large  clumps  of  soil  with  the  probe.  These 
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Dust  Density,  p  gram/ft 


Figure  5.17  Comparison  of  Prediction  with  Experimental 

Observations 


clumps  are  expected  to  exist  close  tc  the  ground.  The  mea¬ 
surement  made  at  the  3-ft  elevation  for  PLUMBBOB  closely 
agrees  in  form  and  magnitude  with  the  pr_ ’icted  results  ex¬ 
cept  for  a  brief  interval.  The  measurement  made  at  the  10-ft 
elevation  is  uniformly  erratic  and  slightly  more  dense  than 
the  corresponding  predicted  value.  The  comparison  for  MET 
shot  is  similar  to  this  latter  comparison.  The  overall 
agreement  is  qualitatively  good.  These  experimental  data  are 
considered  to  be  too  unreliable  to  justify  any  change  in  the 
selected  values  of  the  erosion  constant. 

The  value  of  0.00006  for  the  erosion  coefficient 
was  used  in  the  computations  made  for  the  current  effort. 

A  value  of  0.20  was  used  for  the  transport  constant  K.  These 
values  should  be  sufficient  for  the  current  effort. 

5.4  Typical  Results 

Four  computer  programs  were  written  in  FORTRAN  IV 
for  the  purpose  of  predicting  the  dust  environment  at  the 
elevations  of  interest.  These  programs  are  all  based  upon 
the  previously  developed  cold  dust  models  and  will  also  treat 
the  following  specific  cases: 

(a)  uniform  ground  conditions--single  burst, 

(b)  uniform  ground  conditions --two  or  three  bursts, 

(c)  variable  ground  conditions --s ingle  burst, 

(d)  late -time  dust  -  uniform  ground  conditions  -- 
single  burst. 

These  programs  are  discussed  in  detail  in  Appendix  A.  This 
section  presents  the  results  of  sp-'cific  calculations  and 
illustrates  the  influence  of  some  of  the  pertinent  variables. 
5.4.1.  Blast  Wind  Induced  Dust  Environment 

It  has  been  previously  mentioned  that  the  dust 
density  is  directly  proportional  to  the  value  of  the  erosion 
coefficient.  Furthermore  it  can  be  established  that  some 
surfaces  are  somewhat  nonuniform  and  that  erosion  cannot  be 
expected  to  occur  locally  over  the  entire  surface.  Therefore, 
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a  land -use  factor  has  been  introduced  which  attempts  to 
account  in  a  gross  manner  for  this  character^ tic  of  a  parti¬ 
cular  surface  or  region  of  the  surface.  The  numerical  value 
of  this  factor  lies  between  unity  and  zero  and  must  be  esti¬ 
mated  on  a  highly  subjective  basis.  As  an  example,  an 
area  of  well  sodded  grass  has  been  assigned  a  land-use  factor 
of  0.1  which  implies  that  only  10  percent  of  the  soil  surface 
contributes  to  the  dust  cloud.  In  the  following  results  a 
variety  of  factors  has  been  used. 

Eight  discrete  particle  sizes  are  incorporated  into 
the  computer  programs,  ranging  in  size  from  5  microns  to  5  mm. 
The  dust  density  variation  at  a  site  surrounded  by  a  uniform 
soil  of  a  single  particle  size  was  computed  for  each  of  the 
eight  particle  size  classes.  The  results  are  presented  in 
Figure  5.18.  These  results  apply  at  the  15-ft  level  and  for 
a  relatively  large-yield  weapon.  Figure  5.19  presents  the 
influence  of  the  height  above  the  ground  for  a  typical  site 
in  which  the  soil  distributed  around  the  site  is  not  uniform. 

The  dust  density  will  persist  at  the  site  for  an 
appreciable  period  of  time  after  the  passage  of  the  shock 
wave.  Figure  5.20  illustrates  a  typical  result  for  a  large- 
yield  weapon.  The  positive  phase  duration  of  the  wind  is 
approximately  10  sec.  The  peak  dust  density  occurs  at  ap¬ 
proximately  50  sec.  The  dashed  curve  represents  the  dust 
density  variation  for  discrete  particle  size  groups  and  indi¬ 
cates  when  a  particular  size  class  drops  below  the  height  of 
interest.  The  solid  line  is  merely  a  smooth  curve  illustrat¬ 
ing  the  continuous  decay  that  a  continuous  particle  size 
distribution  would  yield.  The  dust  density  does  not  decay 
below  one-half  the  peak  dust  density  until  approximately 
3  hours  after  the  burst.  It  is  quite  clear  that  the  influence 
of  ambient  winds  will  modify  these  late-time  results,  in 
that  dust  raised  at  other  nearby  regions  will  have  drifted 
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Figure  5.20  Late  Time  Effects 


past  the  site.  For  this  reason  the  late- time  dust  model  was 
developed  and  typical  results  for  late  times  are  presented 
in  the  following  subsection. 

The  variable  ground  dust  code  (single  burst)  per¬ 
mits  one  to  calculate  the  dust  density  variation  at  the  site 
for  the  case  where  the  soil  characteristics  (e.g.,  weight 
factions  and  land  use  factors)  vary  with  ground  range  in 
front  of  and  behind  the  site.  A  series  of  calculations  was 
made  in  which  a  circular  region  R,  around  the  site,  had  a 
stabilized  soil  condition,  such  that  no  dust  could  be  raised 
from  this  region.  The  region  outside  this  circle,  however, 
will  contribute  soil  to  the  dust  cloud.  The  results  of  these 
calculations  are  presented  in  Figure  5.21.  The  maximum  range 
from  the  site,  which  still  contributes  to  the  dust  cloud  at 
the  3ite,  is  approximately  2500  ft. 

The  third  computer  program  permits  one  to  compute 
the  dust  densities  resulting  from  two  or  three  bursts.  These 
buists  must  come  from  the  same  general  direction,  relative  to 
the  site,  and  apply  to  the  uniform  ground  condition  case.  The 
results  for  a  typical  large  yield  case  (for  identical  size 
weapons)  are  illustrated  in  Figure  5.22.  The  additional  dust 
raised  by  the  second  and  third  bursts  contributes  relatively 
little  to  the  dust  density  when  compared  with  the  first  burst. 
This  is  due  to  the  fact  that  although  the  overall  cloud 
height  may  be  nearly  doubled  for  two  bursts,  the  value  of  the 
similarity  variable  is  reduced  by  a  factor  of  approximately 
two.  Such  a  reduction  results  in  a  small  increase  in  the 
dust  density  for  heights  less  than  one-tenth  the  cloud  height 
(see  Figure  5.6). 

5.4.2.  Late -Time  Dust  Environment 

The  overall  dust  cloud  formation  must  first  be  de¬ 
fined  in  general  terms,  and  then  the  apparent  location  of  the 
site  within  the  cloud  at  any  given  time  must  be  established. 

At  this  point  the  blast  wave  parameters  can  be  evaluated  and 
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Figure  5.  21  Influence  of  Soil  Stabilization  Around  Site 
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Figure  5.22  Multiburst  Dust  Cloud 


the  initial  dust  environment  established  at  each  point  of 
interest  in  the  cloud  formation.  It  should  be  noted  that 
for  all  nonsurface  burst  conditions  the  blast  region  will 
be  divided  in  two  shock  reflection  regions;  regular  and  Mach 
reflection  regions.  The  overpressure  variation  with  range 
from  the  origin  for  a  representative  large  yield  burst  con¬ 
dition  is  illustrated  in  Figure  5.23.  The  dividing  point 
between  the  two  reflection  regions  is  defined  as  Rc.  Fig- 
ure  5.24  presents  the  initial  height  of  the  cloud  as  a 
function  of  range  for  the  same  burst  and  yield  conditions. 

The  cloud  height  at  the  origin  vanishes  because  the  horizontal 
air  motion  is  zero.  The  overpressure  is  a  maximum  at  the 
origin  as  indicated  in  Figure  5.24.  The  dust  model  and  the 
associated  blast  wave  idealizations  are  valid  for  low  to 
moderate  overpressure  levels.  At  the  present  time  the  weapons 
effects  subroutines  are  valid  to  overpressures  of  approxi¬ 
mately  60  psi.  For  a  surface  or  near- surface  burst  the  pres¬ 
sure  at  the  origin  will  far  exceed  this  maximum  value  and 
create  a  dust  cloud  environment  which  is  not  well  established. 
This  local  region  cannot  be  handled  properly  at  the  present 
time  and  must  therefore  be  excluded  from  the  present  analysis 
and  calculations.  When  this  region  exists  it  will  be  centered 
about  the  origin  as  indicated  by  the  shaded  area  in  Figure  5.8. 
The  late- time  dust  computer  code  bypasses  this  region  auto¬ 
matically.  It  should  be  noted  that  this  exclusion  will  occur 
only  in  a  very  small  number  of  problems  of  interest. 

Typical  results  are  presented  in  Figures  5.25,  5.26 
and  5.27  for  three  different  surface  wind  magnitudes  and  four 
different  directions  ranging  from  directly  towards  the  origin 
to  directly  away  from  the  origin.  The  site  is  initially 
located  at  the  range  where  the  maximum  cloud  height  exists 
for  these  sets  of  calculations,  hence  the  0=0°  and  9  =  45° 
results  are  characterized  by  cusps.  This  is  due  to  the  pas¬ 
sage  of  the  apparent  path  of  the  site  through  this  cloud 
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Figure  5.27  Influence  of  Wind  Direction,  W  -  10  Knots 


height  maximum  and  therefore  a  local  dust  density  maximum. 

The  zero  dust  density  result  for  0  =  0°  is  due  to  the  passage 
of  the  apparent  site  path  through  the  origin.  The  inter¬ 
mediate  wind  velocity  case  (W  =  4  knots)  appears  to  yield  the 
largest  dust  densities  during  the  time  interval  around  2  to 
3  hours  after  burst.  The  influence  of  location  within  the 
cloud  formation  and  elapsed  time  for  settling  combine  to  yield 
this  approximate  ’’worst"  case  situation. 

The  influence  of  the  wind  velocity  is  demonstrated 
in  Figure  5.28  for  a  typical  burst  situation  and  wind  direc¬ 
tion.  The  worst  case  during  the  first  four  hours  is  clearly 
given  by  a  wind  velocity  of  0.5  mph.  The  larger  surface  winds 
clearly  blow  the  dust  away.  At  latter  times  the  worst  case 
corresponds  to  a  wind  velocity  of  0.02  mph.  This  case  yields 
greater  dust  densities  than  for  zero  wind  velocity  because 
the  wind  magnitude  is  just  great  enough  to  keep  airborne 
(C  <  0)  the  smallest  particle  size  class  used  (5  microns) 
Thus  there  is  a  real,  although  somewhat  subtle  influence,  of 
particle  size. 

The  influence  of  weapon  size  or  yield  is  presented 
in  Figure  5.29.  Here,  the  effect  is  rather  clearcut  and 
significant.  Figure  5.30  illustrates  the  influence  of  height 
of  burst  on  the  dust  density.  The  only  large  difference  is 
for  the  near  overhead  case  of  a  height  of  burst  of  30,000  ft, 
and  the  difference  occurs  only  during  the  first  four  hour 
period.  Figure  5.31  presents  the  influence  of  elevation 
above  the  ground  surface.  Clearly  very  little  difference 
in  dust  density  exists  in  the  region  of  from  10  to  30  feet 
above  the  ground  surface.  Figure  5.32  presents  the  results 
of  varying  the  blast  overpressure  at  the  site.  The  late 'time 
dust  environment  is  not  very  sensitive  to  this  blast  variable. 


177 


Yield  -  25.0  MT 


5.29  Influence  of  Yield 
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Figure  5,32  Influence  of  Overpressure 


5.5 


List  of  Symbols 
a 


B 

c 

C 

Cd 

D 

f 

k 

K 

M 

N 

n 


P 

Re 

R, Rq , Ri ,R2 , . . R4  >  RS s  Rt 
Rc 

t » tf, • . . t5 

ti 

to 


Erosion  coefficient 

Constant 

Constant 

Constant 

Constant 

Drag  coefficient 

Particle  diameter 

Factor 

Cloud  height 

Exponential  coefficient 

Transport  coefficient 

Weight  of  dust 

Number  of  particle  diameters 

Number  of  particles 

Weight  Fraction 
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Density  of  particles 

Specific  density  of  soil 
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5.6  Dust  Cloud  Generated  by  Crater  Ejecta 

Several  sources  contribute  to  the  generation  of  a 
dust  cloud  after  the  detonation  of  a  nuclear  weapon.  One  such 
source  is  crater  ejecta.  An  estimate  of  the  quantity  of  dust 
generated  via  crater  ejecta  must  be  made  and  compared  to  the 
contribution  of  other  sources.  A  previous  attempt'*' "*  to  esti¬ 
mate  crater  ejecta  dust  assumed  that  all  the  dust  was  distri¬ 
buted  over  a  cloud  height  of  25  ft.  This  assumption  is 
arbitrary  and  indicated  that  the  resulting  dust  cloud  density 
was  very  high.  The  procedure  used  to  make  these  estimates 
did  not  allow  the  crater  ejecta  to  reach  the  location  of 
interest  over  a  finite  period  of  time. 

The  characteristics  of  a  dust  cloud  generated  by 
crater  ejecta  are  very  complex.  However,  a  simple  model  that 
characterizes  the  phenomenon  should  be  adequate  to  predict 
the  dust  cloud  density  provided  the  predictions  are  conserva¬ 
tive  and  are  small  in  comparison  to  the  other  dust  cloud 
sources . 

C  fi 

It  has  been  established"*’  that  the  crater  ejecta, 
which  are  thrown  distances  greater  than  several  crater 


184 


diameters,  originate  at  or  near  the  soil  surface.  It  is 
transported  ballistically  as  a  result  of  ito  initial  velocity. 
The  air  is  also  being  accelerated  at  this  time.  However, 
at  late  times  the  high  velocity  winds  outrun  the  crater 
ejecta  particles.  A  model  which  would  adequately  describe 
this  phenomenon  is  one  in  which  (1)  the  quantity  and  initial 
flight  conditions  of  the  ballistic  ejecta  are  given,  and 
(2)  the  time  \istory  of  the  ejecta  particles  are  followed 
through  the  blast-accelerated  air  until  they  settle  at  the 
point  of  interest.  The  evaluation  of  such  a  mode)  would 
require  considerable  effort  if  the  influences  of  the  blast 
winds  and  drag  forces  are  included.  It  would  be  impossible 
to  determine  which  particles  reach  a  given  range  without 
examining  the  trajectories  of  all  of  the  particles  of  the 
crater  ejecta.  Therefore,  for  the  purpose  of  making  an  ini¬ 
tial  estimate  of  the.  dust  density,  the  influence  of  both  air 
motion  and  the  drag  forces  acting  on  the  crater  ejecta  has 
been  neglected.  The  time  of  flight  t£  for  a  specified 
quantity  of  debris  to  reach  the  range  distance  R  is  given  by 

2V  Tl  /»  \“ 

tf  =  “g2  sin  7  arc  sin  -|j  ,  (5.52) 

\Vo  i  _ 

where 

t £  =  time  of  particle  flight, 

R  =  range , 

VQ  =  initial  particle  velocity,  and 
g  =  gravitational  constant. 

The  initial  angle  of  flight  is  uniquely  related  to  the  range 
and  initial  velocity  of  the  particles,  such  that  one  pair  of 
velocity-initial  angle  values  will  carry  the  particle  to 
the  range  of  interest.  This  angle  0  is  given  as 
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(5.53) 


9 


arc  sin 


It  should  be  noted  that  9  can  be  measured  either  from  the 
horizontal  or  vertical  axis  since  the  relationship  is  sym¬ 
metric  about  9  =  tt/4.  Furthermore,  since  the  value  of  the 
arc  sin  is  bounded,  there  exists  a  minimum  value  of  the 
initial  velocity  which  will  result  in  a  flight  of  range  R. 

This  value  occurs  at  9  =  rr/4.  Experimental  evidence  from 
crater  ejecta  studies'****  indicates  that  the  maximum  initial 
velocity  of  the  crater  ejecta  is  approximately  1400  fps.  Thus, 
it  is  possible  to  place  limits  upon  the  value  of  the  initial 
velocity  of  flight.  The  trajectories  of  the  crater  ejecta 
are  illustrated  in  Figure  5.33. 


The  quantity  of  crater  ejecta  that  reaches  a  given 
range  is  subject  to  some  uncertainty.  The  results  of 
Carlson"* *^  have  been  used  for  these  calculations.  The  weight 

of  ejecta  w,  per  unit  area,  that  reaches  a  specified  range 
is  given  by 

IK  \3*86 

w  =  -i  D_  i  «-~l  .  (5.54) 


„  ! 


1 


where 


w  =  weight  of  crater  ejecta  per  unit  area, 
Da  =  apparent  crater  depth, 

Ra  =  apparent  crater  radius,  and 
Ps  =  particle  density. 


The  manner  in  which  the  crater  dimensions  scale 
with  weapon  yield  W  is  questionable.  The  weight  of  crater 
ejecta  scales  as  or  W~"^.  The  estimates  made  herein 

use  the  latter  scaling  law. 
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The  distribution  of  ejecta  mass  with  initial 
velocity  is  by  no  means  clear.  It  would  appear  reasonable 
to  assume  that  the  majority  of  the  mass  leaves  the  crater 
with  a  low  velocity  and  that  very  little  ejecta  mass  has  an 
initial  velocity  of  1400  fps.  This  computation  assumes  a 
linear  relationship  between  mass  and  initial  velocity,  the 
peak  occurring  at  the  minimum  velocity,  and  a  value  of  zero 
occurring  at  VQ  =  1400  fps.  This  mass  is  distributed  symmet¬ 
rically  about  9  =  n/4.  The  corresponding  polar  diagram  is 
illustrated  in  Figure  5.34  for  a  large  weapon  yield.  Based 
on  the  above,  it  is  possible  to  compute  the  weight  of  ejecta 
(per  unit  area)  that  reaches  the  range  of  interest  in  a 
particular  time  interval. 

The  density  of  the  dust  due  to  the  crater  ejecta 
can  be  obtained  by  dividing  the  time  rate  of  mass  (weight) 
deposition  by  the  vertical  velocity  of  the  ejecta.  This 
velocity  is  VQ  cos  0 .  It  should  be  noted  that  the  ejecta  is 

impacting  on  the  ground  at  its  initial  velocity,  since  no 

drag  effects  were  incorporated  into  the  model.  One  would 

expect  the  ejecta  to  fall  at  a  lower  velocity,  such  as  its 

associated  terminal  velocity.  The  results  of  calculations 

based  upon  the  above  model  are  presented  in  Figure  5.35. 

They  are  given  as  a  function  of  time  for  a  small  and  large 

weapon  yield.  The  peak  dust-cloud  density  for  the  smaller 

yield  is  approximately  0.66  gr/ft^  and  is  considerably  less 

than  the  dust-cloud  densities  predicted  for  the  blast-induced 

dust  clouds.  The  peak  dust-cloud  density  for  the  large  yield 

3 

weapon  was  less  than  0.2  gr/ft  . 


Figure  5.35  Dust  Density  due  to  Crater  Ejecta 


It  is  difficult  to  ascertain  that  tl~e  above  calcu¬ 
lations  are  indeed  conservative,  although  it  is  obvious  that 
the  actual  times  of  flight  will  be  considerably  greater  than 
those  computed  by  the  previous  model.  The  above  model  indi¬ 
cates  that  the  debris  arrives  at  the  site  between  19  and 
85  seconds  after  the  weapon  detonates. 
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EXTRACTED 
Chapter  6  (U) 

AIR-TEMPERATURE  STUDIES 

6.1  Introduction 

The  thermal  radiation  emitted  by  the  fireball  of 
an  atmospheric  burst  will  be  transmitted  through  the  sur¬ 
rounding  air  with  relatively  little  absorption  until  it 
reaches  some  relatively  opaque  substance  such  as  the  ground, 
or  perhaps  a  dust  cloud.  The  deposition  of  the  thermal 
v energy  in  the  opaque  material  will  raise  its  temperature 
locally. 

In  the  absence  of  a  dust  cloud,  some  of  the  thermal 
energy  will  heat  the  ground  surface  (and  some  portion  immedi¬ 
ately  below  the  surface)  to  relatively  high  temperatures. 

The  soil  may  respond  mechanically  to  this  thermal  pulse, 
depending  upon  the  temperature  and  thermal  gradients  which 
are  generated.  The  soil  may  also  melt  or  vaporize  in  some 
instances . 

The  passage  of  the  shock  wave  at  a  particular  posi¬ 
tion  on  the  ground  (say  the  site  location)  will  result  in  the 
start  of  the  blast-induced  winds  at  that  position.  If  the 
soil  is  such  that  it  can  be  eroded  by  the  blast  winds,  a 
dust  cloud  will  begin  to  form.  This  phenomenon  was  described 
in  the  preceding  chapter  for  the  case  of  cold  dust  (i.e.,  no 
thermal  effects  were  considered) .  As  a  result  of  the  deposi¬ 
tion  of  the  thermal  energy  in  the  soil,  some  of  the  particu¬ 
late  matter  which  forms  the  dust  cloud  will  be  quite  hot 
(relative  to  the  air) .  The  air  will  be  heated  by  the  shock 
compression  and  then  cooled  as  the  gas  begins  to  expand. 

The  passage  of  the  hot  dust  through  the  air  will  result  in 
the  exchange  of  thermal  energy  between  the  dust  and  the  air. 
Furthermore,  additional  energy  will  be  deposited  in  the  dust 
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cloud  by  the  process  of  absorption  of  the  fireball  thermal 
radiation.  This  energy  will  also  be  transmitted  to  the  air, 
thus  further  altering  the  air  temperature  at  the  heights  of 
interest . 

The  overall  air-temperature  problem  is  quite  complex 
and  considerable  simplification  must  be  resorted  to  in  order 
to  achieve  an  engineering  estimate  of  the  air-temperature 
variation  at  the  potential  site  locations.  The  air-tempera¬ 
ture  model  which  is  described  in  this  chapter  is  based  upon 
the  dust  model  described  in  the  preceding  chapter.  All  of 
the  assumptions  and  simplifications  made  for  the  dust  model 
are  used  herein,  and  the  basic  assumption  relating  to  the 
uncoupling  of  the  gas  dynamic  field  from  the  local  dust  and 
air  temperature  variations  is  retained. 

The  following  section  presents  a  description  of  the 
thermal  source  as  it  relates  to  the  current  problem.  Sec¬ 
tion  6.3  deals  with  the  deposition  of  the  thermal  radiation 
in  the  soil  and  the  resulting  transient  soil  temperature. 
Section  6.4  treats  the  deposition  of  the  thermal  energy  in 
a  dust  cloud.  The  basic  IITRI  air-temperature  model  which 
incorporates  all  of  the  preceding  effects  together  with  the 
dust  cloud  characteristics  is  then  presented  in  Section  6.5. 
The  next  two  sections  (Sections  6.6  and  6.7)  present  typical 
results  for  single-and  multiple -burst  situations,  respective¬ 
ly,  and  illustrate  the  influence  of  certain  variables.  The 
last  section  presents  a  list  of  symbols  used  in  this  chapter. 

6.2  Thermal  Radiation  from  the  Fireball 

The  thermal  energy  emitted  by  the  fireball  is 
6  1 

reasonably  well  known.  *  For  the  purpose  of  this  study,  a 
simple  representation  of  the  thermal  emission  will  be  used. 

The  fireball  will  be  treated  as  a  point  source  rather  than 
the  finite  source  it  actually  is.  This  approximation  is 
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reasonable  for  distances  much  greater  than  the  fireball 
diameter.  This  is  generally  the  case  for  the  current  problem. 


The  fireball  emission  is  treated  as  black-body 
radiation  at  an  effective  temperature  T  .  The  thermal  flux 
Fsis 

Fs=oTe4  (6.1) 

where  rj  is  the  Stefan's  constant.  The  total  radiant  power 
emitted  from  a  fireball  of  radius  R  is 


P 


4"R2aT  4  . 
e 


(6.2) 


The  flux  at  a  slant-range  distance  D  from  the  fireball 
(R«D)  is 

F  =  t  P/4ttD2  (6.3) 

a 

where  t  is  the  atmospheric  transmittance.  The  radiant  flux 
intensity  ttB^ ,  as  a  function  of  the  fireball  surface  tempers 
ture  T ,  is  shown  in  F igure  6.1. 

For  a  weapon  yield  of  W  kilotons,  the  approximate 

6  1 

maximum  thermal  emission  ’  (calories  per  sec)  is 

Pmax  ~  4(1()12>  Vw  .  (6.4) 

The  approximate  time  (sec)  after  explosion  for  the  tempera¬ 
ture  maximum  is 

‘max  ~  °-032  Vw"  •  (6.5) 


An  alternate  expression  involving  the  height  of  burst  is 


‘max  =  °-82  W42  exp (-.09  H„)  . 


(6.6) 
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Radiant  Power,  ttB,  >  Ergs/cm  -sec 


Wavelength,  X,  Angstroms 


Figure  6.1 


Radiant  Power  of  a  Black  Body  as  a  Function  of 
Wavelength  at  Various  Temperatures 
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This  maximum  is  for  the  second  thermal  pulse;  the  first  pulse 
is  of  such  brief  duration  that  its  thermal  contribution  is 
neglected. 

The  approximate  scaled  emission  as  a  function  of 
scaled  time  is  shown  in  Figure  6.2. 

Some  absorption  occurs  in  the  atmosphere,  so  the 
radiant  intensity  is  attenuated  by  the  factor 

Ta  =  exp(-KD),  (6,7) 

where  i  is  the  atmospheric  transmittance  and  k  is  an  ab¬ 
sorption  coefficient  averaged  over  the  wavelength  spectrum 

and  the  slant  range  D.  Typical  values  of  transmittance  as  a 

6  1 

function  of  slant  range  are  shown  in  Figure  6.3  (Glasstone  *  ). 
The  height  of  burst  is  Hg  =  D  cos  9,  so  that 

t  =  exp(- k  Hg  sec  9 )  .  (6.8) 

6  9 

Gibbons  *  tabulates  the  values  of  k Hg  (which  he  calls  the 
optical  thickness)  as  a  function  of  Hg  for  radiation  of  wave¬ 
length  0.65  micron. 

6.3  Soil  Temperature 

6.3.1  General 

The  thermal  energy  which  is  incident  upon  the  soil 
surface  will  be  absorbed  at  the  surface  and  will  propagate 
downward  by  the  process  of  thermal  conduction.  Furthermore t 
if  the  surface  temperature  becomes  quite  hot  it  will  re-radi¬ 
ate  energy  back  up  into  space.  This  energy  deposition  will 
manifest  itself  as  a  temperature  variation  in  the  soil  near 
the  surface  of  the  soil.  When  the  shock  wave  sweeps  past 
the  soil,  the  soil  will  be  eroded  by  the  blast  winds,  and 
soil  particles  of  various  temperatures  will  enter  the  dust 
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cloud  at  various  times.  Thus,  it  is  necessary  to  determine 
the  vertical  temperature  distribution  in  the  soil  at  the 
time  of  arrival  of  the  shock  wave  at  the  site. 

The  deposition  of  the  thermal  energy  in  the  soil 
represents  a  mechanism  for  concentrating,  locally,  a  signifi¬ 
cant  portion  of  the  energy  (per  unit  area)  released  by  the 
weapon.  This  concentration  of  energy  can  be  maximized  by 
neglecting  any  mechanisms  which  might  tend  to  distribute 
the  energy  more  uniformly.  Therefore,  any  mechanical  response 
that  the  soil  might  experience  will  be  neglected.  Also,  any 
changes  of  phase  of  the  soil  or  its  constituents  (such  as 
water)  will  be  neglected. 

The  soil  temperature  distribution  at  the  time  of 
arrival  of  the  shock  wave  is  determined  by  the  solution  of 
a  nonsteady  one -dimensional  heat  conduction  problem  subject 
to  the  thermal  flux  from  the  weapon,  including  any  re-radi¬ 
ation  which  might  occur. 

6.3.2  Conduction  Solution 

The  thermal  energy  release  of  a  nuclear  weapon 
creates  a  temperature  distribution  T(x,t)  in  the  soil  ac¬ 
cording  to 


and 


where 


T(x,0)  -  TqJ 


\ 


dT 

Sx 


x=0 


=  f(t)  =  -^(t)  -  Qr (t) 


(6.9) 


x  =  distance  below  the  surface  of  the  soil, 
t  =  time, 

k  =  thermal  conductivity  of  the  soil,  and 
&  =  thermal  diffusivity  of  the  soil. 
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The  incident  thermal  power  density  is  / 

in  t a  cos  ® 

Q.  (t)  =  3. 994467 (10iU)  -2 - , -  P  (t)  ,  (6.10) 

i  Dz 

and  the  re-radiated  thermal  power  density  is 

Qr(t)  =  a  T4(0,t)  -  Tq4  (6.11) 

where  a  =  Stefan's  constant.  The  quantities  t  ,  W,  9,  and 

a 

D  are  computed  as  weapons  effects  variables;  Pe(t)  is  the 
scaled  fireball  thermal  power  provided  in  tabular  form  as 
a  function  of  time;  and  T(0,t)  is  the  time  varying  (absolute) 
surface  temperature  which  is,  itself,  the  solution  to  Equa¬ 
tion  6.9.  t  was  defined  previously  as  the  atmospheric 
attenuation  factor,  W  is  the  weapon  yield  in  kilotons,  9  is 
the  angle  of  incidence  of  the  radiation  upon  the  soil  surface, 
and  D  is  the  slant  range  from  the  burst  point  of  the  weapon. 

Tq  is  the  initial  temperature  of  the  soil. 

The  general  solution  to  Equation  6.11  is 


.(6.12) 


Letting  g(T  )  =  -  f(T)  =  CjP^t)  -  C2  T4(0,t)  -  Tq4  ,  Equa¬ 

tion  6.12  is,  at  the  surface  (x=0) , 


AT (0,t) 


and  below  the  surface  (x  >  0) , 


(6.13) 


(6.14) 
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An  iterative  solution  of  Equation  6.13  is  necessary,  since 
g(T)  is  a  function  of  AT(0,t),  but  once  its  solution  is  ob¬ 
tained,  Equation  6.14  can  be  solved  explicitly  for  AT(x,t) 
(i.e.,  at  any  depth  x)  without  iteration. 

Since  the  integrand  of  Equation  6.13  is  unbounded 
as  T—^t,  and  the  integrand  of  Equation  6.14,  although 
everywhere  bounded,  has  a  sharp  spike  for  values  of  t  near 
t  and  small  x,  standard  numerical  integration  schemes  would 
require  considerable  refinement  for  accurate  solutions. 

Thus,  an  analytic  method  of  solution  was  devised  for  each 
of  these  equations . 

Figure  6.4  is  a  graph  of  g(t)  where  it  is  assumed 
that  ATCOjt^)  and  therefore,  g^  =  g(t^)  is  known  at  all  t^ 
for  i<m.  The  g(-r)  curve  can  be  approximated  by  straight 
line  segments  over  suitably  small  increments  of  time,  i.e., 

g(T )  =  t  +  bi,  ti^T  <  ti+1  ,  (6.15) 


Substituting  Equation  6.15,  the  integration  of  Equation  6.13 
can  be  approximated  by 
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g(t) 


AT (0,t  )  =  l\am  A2t  +  t  ,)  +  3b  n  ]  Vt  -  t 

m'  JL  m-1'  m  m-17  m-lJ  f  m  m-1 

m-2 

+  Y"  4  (2a.  t  +  a.t.  +  3b.)  Vt  -  t. 
J  s  l  m  ii  i7  *  m  i 

i=l 

-  (2a. t  +  a.t.,!  +  3b.)  V"t  -  t.,..  1 

i  m  l  i+l  i7  *  m  l+l  -» 


(6.17) 


where  a  ,  and  b  ,  are  functions  of  T(0,tm).  Assuming  the 
solution  has  previously  been  obtained  at  times  tm  prior  to 
the  current  value  of  tm,  the  summation  above  (to  be  denoted 
below  by  S)  is  known.  Substituting  Equations  6.16  and  the 
definition  of  g(tm)  into  Equation  6.17  gives  the  fourth 
order  algebraic  equation 

T4(0,t:m)  +  C3  T(0,tm)  -  C4  =  0  (6.18) 


where 


C3  =  1/^  C2  -V  tm  -  VP  and 

4  Sm-l  .  ClPe(tm) 

C4  =  (T0  +  S)  C3  +  Tq  +  ^  +  -  C—  • 

The  solution  to  Equation  6.18  can  be  found  from  tie  simple 
recursion  formula 


(n+1)  (0,tm)  = 


3[T(n)(0,tm)]+C, 


4fT(n)  (0.t_)1 


(6.19) 


+  r.. 


where  the  superscript  n  represents  the  n£h  iteration,  and 
initially,  T(1)(0,tm)  is  taken  as  TKO,^). 
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Equation  6.19  is  thus  solved  at  each  tffi  successive¬ 
ly,  as  time  increases,  resulting  in  both  T(0,t)  and  g(t) 
or,  equivalently,  the  piece-wise  linear  parameters  a^  and  b^. 

The  subsurface  temperatures  can  now  be  computed 
directly  by  integration  of  Equation  6.14  according  to  the 
approximation 


AT(x,tm) 


where  the  ai  and  b^  are  known.  Thus, 


for  the  soil  temperature  at  any  depth  x  and  time  t^. 

6.3.3  Temperature  Distribution  in  Soils 

A  computer  subprogram  was  written  in  FORTRAN  IV 
in  accordance  with  the  previous  analysis  to  determine  the 
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temperature  distribution  in  the  soil  at  the  time  of  arrival 
of  the  shock  wave. 

Since  soils  are  both  varied  and  complex,  it  was 
necessary  to  investigate  the  pertinent  thermal  properties  for 
representative  soils  and  to  investigate  the  influence  that 
any  uncertainty  of  these  properties  has  on  both  the  tempera-  * 

ture  distributions  and  the  variation  which  might  be  intro¬ 
duced  into  the  air-temperature  prediction. 

Table  6.1  presents  a  brief  summary  of  the  thermal 
6  3 

properties  *  of  soil.  The  variation  is  not  too  large  except 
for  a  few  extreme  cases  such  as  mud  or  very  wet  soils.  A 
series  of  temperature  distributions  were  computed  for  a 
number  of  cases  for  both  a  small  yield  and  a  large  yield 
case.  These  temperature  distributions  are  presented  in 
Figures  6.5  and  6.6.  The  temperature  distributions  for  the 
large  yield  weapon  were  used  in  a  series  of  air -tempera ture 
calculations.  These  results  will  be  discussed  in  a  later 
section  of  this  chapter. 

A  low  value  for  the  thermal  conductivity  yields, 

as  is  expected,  higher  surface  temperatures.  The  total 

energy  deposited  in  the  soil  increases  with  both  a  decreasing 

value  of  the  thermal  conductivity  and  an  increasing  value 

of  thermal  dif fusivity .  In  the  absence  of  specific  site  data, 

a  value  for  the  thermal  conductivity  of  0.10  Btu/hr-ft-°F 

is  recommended  for  general  use.  Similarly,  a  value  of 
o 

0.007  ft  /hr  r  the  thermal  diffusivity  is  recommended. 

6.4  Radiation  Absorption  in  Dust  Cloud 

The  radiation  of  energy  from  the  above-ground 
fireball  irradiates  the  dust  cloud  generated  at  the  ground 
surface.  The  amount  of  energy  radiated  is  a  function  of 
wavelength  and  fireba.i.  temperature.  For  instance,  it  is  we  11 
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Table  6.1  THERMAL  PROPERTIES  OF  SOIL 


Temperature 


120Q 


Figure 


Depth,  z,  inches 

6-5  Soil  Temperature  Distribution,  Small  Yield 
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known  that  most  of  the  energy  is  radiated  in  the  visible 
range  of  0.3  to  0.7  micron  wavelengths  for  a  5000°K  fireball. 

In  the  infrared  rarge  of  0.7  to  15  micron  wavelengths,  and 
in  the  X-ray  range  of  0.001  to  0.015  micron,  there  is  a 
two-order  magnitude  of  radiant  energy  reduction. 

The  transmission  of  radiant  energy  through  a  cloud 
of  dust  particles  which  are  distributed  randomly  in  space 
and  time,  and  which  have  a  uniform  diameter  (monodisperse) , 
can  be  calculated  from  the  formula 

I  =  Io  exp(-i  X)  (6.21) 

where  I  is  the  incident  radiation  flux  at  the  top  of  the 
o 

dust  cloud,  I  is  tha  radiant  flux  in  the  dust  cloud  at  a 
distance  "X"  from  the  top  of  the  dust  cloud,  and  a  is  the 
absorption  coefficient. 

The  absorption  coefficient  can  be  calculated  on 
the  basis  of  the  dust  cloud  being  composed  of  spherical 
particles,  with  an  average  separation  between  spherical 
particles  being  much  greater  than  the  wavelength  of  the 
incident  radiant  energy.  Thus,  the  scattering  and  absorption 
properties  of  the  cloud  can  be  determined  in  terms  of  the 
scattering  and  absorption  cross  sections  of  a  single  spheri¬ 
cal  particle. 

The  absorption  and  scattering  of  a  single  spherical 
particle  have  been  evaluated  analytically  by  Mie.  A  plane 
electromagnetic  wave,  incident  upon  a  spherical  particle, 
excites  the  various  natural  electric  and  magnetic  resonance 
modes  in  -the  particle.  Thereafter,  the  particle  re-radiates 
a  portion  of  its  energy  in  all  directions  and  absorbs  the 
other  portion  of  the  energy,  at  the  same  frequency  as  the 
incident  energy.  Mie’s  equations  have  been  solved  numerically, 
and  the  computer  results  are  often  presented  in  terms  of  a 
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normalized  cross  section  Q,  referring  to  the  spherical  particle 

2 

cross-sectional  area  ^r  .  These  cross  sections  are  a  function 
of  the  size  parameter  2tti:/\  where  r  is  the  radius  of  the 
spherical  particle,  and  X  is  the  wavelength  of  the  incident 
energy  for  each  value  m  of  complex  index  of  refraction  of 
the  particle  material. 

For  very  large  values  of  2nr /X ,  in  the  range  of 
2rrr/X  >  50,  the  sum  of  the  normalized  absorption  (Q0)  and 
scattering  (Q  )  cress  sections  approaches  the  value  2. 

For  very  small  values  of  2rrr/X  ,  in  the  range  of 
2nr/X  <  0.1,  the  scattering  and  absorption  cross  sections  are 
related  to  2rrr/x  ,  raised  to  an  integral  power. 

For  the  range  of  values  of  2rrr/x  between  4  and  50, 
there  is  a  strong  oscillatory  relation  between  Q  and  2nr/a. 

It  has  become  customary  to  refer  to  the  normalized 
cross  sections  Q  as  efficiencies.^*^  Also,  the  sum  of  the 
absorption  and  scattering  efficiencies  is  called  the  extinc¬ 
tion  efficiency. 


‘abs 


+  Q  =  Q  . 

xsca  Xext 


Typical  results  '  of  those  efficiency  (Q)-versus-size  param¬ 
eter  (2nr/X)  curves  for  a  particular  complex  value  of  index 
of  refraction  (m)  are  shown  in  Figures  6.7  through  6.10. 

2 

Consider  that  I  is  the  intensity  (watts /u  )  of 
the  incident  radiation  which  is  not  polarized.  A  spherical 
particle  will  intercept  Qextnr  IQ  (watts)  from  the  incident 
beam,  where  Qabswr^0  will  be  the  absorbed  energy  and 
Q  nr^I  will  be  the  scattered  energy  of  the  particle. 

The  scattered  energy  has  preferred  directions 
depending  upon  the  size  parameter  and  index  of  refraction, 
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Figure  6.7  Extinction  Efficiency  for 
Greater  than  1.0 :m  ■  1.50 
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Figure  6.9  Extinction,  Scattering,  and  Absorption 

Efficiencies  vs  Size  Parameter:  m  =  1,315 
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and  has  a  polarization  effect.  For  instance,  when  2rrr/X<<l, 
the  well-known  Rayleigh  Law  for  the  scattering  of  unpolar¬ 
ized  radiation  is 


where  I  is  the  intensity  of  the  scattered  radiation  in  the 
direction  0,  V  is  the  volume  of  the  particle,  and  R  is  the 
distance  between  the  particle  and  the  observer.  The  radi¬ 
ation  intensity  scattered  by  the  particle  consists  of  two 
incoherent  plane-polarized  components,  i^  and  whose 
planes  are  mutually  perpendicular.  The  factor  of  unity  in 

o 

the  first  parenthesis  of  IQ  governs  i^;  while  the  factor  cos  0 
in  the  last  parenthesis  of  IQ  governs  i ^ •  At  0  =  rr/2,  the 
cos^0  factor  is  zero  and  the  scattered  radiation  becomes 
completely  plane-polarized,  in  a  plane  perpendicular  to  the 
plane  of  observation. 

The  angular  distribution  of  intensity  is  symmetri¬ 
cal  about  a  plane  normal  to  the  incident  light,  i.e.,  as 
much  light  is  scattered  forward  as  backward  when  2rrr/X<0.1. 
However,  about  a  1000-fold  or  more  ratio  exists  favoring 
the  forward  over  the  backward  scattering  for  2TTr/X>l.  Since 
the  average  size  dust  particle  is  12.5  microns  in  radius, 
and  we  are  concerned  with  visible  light  in  the  range  of 
0.5  micron,  2TTrA>l  and  there  is  strong  forward  scattering. 
The  sensitivity  of  the  scattering  to  the  index  of  refraction 
m  and  size  parameter  a  is  shown  in  the  following  diagram. 

Note  that  the  total  scattered  energy  is  given  by  i^  + 

The  particles  in  a  dust  cloud  scatter  energy  mostly 
in  the  forward  direction,  and  weakly  to  the  backward  and 
sideward  directions  as  shown  in  Figures  6.11  and  6.12.  The 
backward  scattered  energy  will  be  small  and  appear  as  re¬ 
flected  energy  from  the  dust  cloud  at  the  top,  and  the 
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Figure  6.11  Polar  Diagram  of  S 
Components  for  0,1 
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Incident  Light 
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Figure  6.12  Polar  Diagram  of  Scattering  for  i^  and  i ^ 
Components  .or  C.33p  Radius  Oil  Droplet 
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forward  scattered  energy  will  be  large  and  appear  as  trans¬ 
mitted  energy  in  the  interior  of  the  dust  cloud.  This  is 
the  absorbed  energy  of  the  particles  that  must  be  used  to 
calculate  the  temperature  changes  in  the  dust  cloud. 

6  7 

Kattawar  *  calculates  Qabs>  for  several  materials 
m,  for  size  parameters  2TTr/X  =  0.1,  1.0,  and  10.  These  charts 
can  be  used  to  determine  the  energy  absorbed  at  various  dis¬ 
tances  into  the  cloud  (X)  according  to: 

I  =  IQ  exp(-NQabs  t  r2X) 

where  N  is  the  number  of  particles  per  cubic  centimeter, 

Qabs  is  the  absorption  efficiency  factor,  and  r  is  the  radius 
of  the  spherical  particle.  Since  the  number  of  particles 
per  cubic  centimeter  depends  upon  the  density  of  the  dust 
cloud  p,  the  density  of  the  dust  particle  material  p  ,  and 
the  radius  of  the  particle  r,  according  to 


N  = 


Dc  1 


the  intensity  of  the  energy  in  the  cloud  at  a  particular 
depth  is 

T  T  3  n  P  x\ 

1  =  *0  exP  -  %  Qabs  p-  r  * 

'  c  I  ^ 

The  absorption  efficiency  depends  upon  the  wave¬ 
length  of  the  incident  radiation  X  ,  via  the  size  parameter 
2nr/X .  The  wavelengths  of  interest  are 

X  *  0.280  —>  0.397  i-i  (ultraviolet), 

X  =  0.397  — >  0.700^1  (visible),  and 

X  =  0.700  —5*  15  — >  150  u  (infrared). 
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Assuming  that  the  average  size  particle  is 
r  =  12.5  (microns),  then 


2rrr 

X 


=  ^  =0.1  (electromagnetic). 


2rrr  _  2 

t — r 

2nr  _  2 

t - r 


=i.o  (infrared), 

=  10.0  (infrared). 


T1  =  (0^785u)^  =  100 -°  (visible  red),  and 

=  |.0C12j5u)  =  200.0  (visible  violet). 

In  the  visible  region  where  the  fireball  emits 
most  of  its  energy. 


2^r  >100 

A. 

and 

Qabs  V  ^  • 


Thus,  the  energy  absorption  in  the  dust  cloud  simplifies  to 


I 


exp 


3 

In  the  special  case  where  the  dust  cloud  is  P  =  0.01  lb/ft  , 
and  each  dust  particle  has  a  specific  gravity  of  2.65  and 
a  radius  of  12. 5|i , 


I  =  IQ  exp 


X 

TO 


where  X  is  measured  in  inches .  Thus ,  the  intensity  of  the 
radiation  is  reduced  by  approximately  2.7  every  10  inches 
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t 


\ 


3 

(25.4  cm).  Since  there  are  7,344  particles/cm  ,  there  is  an 
average  spacing  of  0.051  cm  between  particles.  The  optical 
depth  is  increased  beyond  10  in.  when  o  is  increased,  p  is 
decreased,  and  r  is  increased. 

In  the  infrared-to-electromagnetic  region  where  the 
fireball  emits  a  small  fraction  of  its  energy, 

0.1<~£  <  10.0 

A 

and 

Qabs  <  1*°  * 


Furthermore,  Q  ^  becomes  a  function  of  the  size  parameter 
2nr/\  and  the  complex  index  of  refraction  m,  as  shown  in 
Figures  6.13  through  6.15.  Consider  the  following  table. 

Table  6.2  ABSORPTION 


L  _ 1 

Index  of 
Refraction 

^absorption 

Type 

Size  Parameter 

m=n^-in2 

a  =  10 
Infrared 

a  =  1.0 
for  r=25u 

a  =  0.1 

E.M.  for  r=25u 

1 

1.33-1. OOOi 

0.4500 

1.5000 

0.2100 

2 

1. 33-0.100i 

0.8000 

0.2800 

0.0210 

3 

1.33-0«.010i 

0.3500 

0.2800 

0.0021 

4 

1.33-0.001i 

0.0410 

0.0028 

C.0002 

5 

10-10. OOOi 

0.3000 

0.4500 

0.03000 

6 

10-1. OOOi 

0.4500 

0.8200 

0.00550 

7 

10-0. 1001 

0.4500 

0.2000 

0.00055 

8 

10-0. OlOi 

0.3000 

C. 0200 

0.0006 

9 

10-0. OOli 

_ 

0.0500 

0.0020 

0.00000 

: .  ... - - ; 
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Figure  6.15  Absorption  Coefficient 
Index  of  Refraction,  2 
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For  types  1,  2,  5,  and  6  dust  materials,  the  previ¬ 
ous  numerical  example  which  resulted  in  a  10-in.  e -folding 
distance  for  radiation  in  the  visible  region  is  also  appli¬ 
cable  to  well  within  an  order  of  magnitude  for  the  visible 
infrared  region,  but  must  be  strongly  increased  for  the 
inf rared-to-e lectromagnet ic  region . 

Material  1  corresponds  to  water  in  the  optical 
region,  while  material  6  corresponds  to  water  in  the  elec¬ 
tromagnetic  region.  Material  2  corresponds  to  metals  in 
the  optical  region,  while  material  5  corresponds  to  metals 
in  the  infrared  region. 

Consider  a  fireball  temperature  of  5000°K  which 
results  in  most  of  the  emitted  energy  being  in  the  visible 
range  of  radiation,  and  an  average  dust  size  particle  of 
r  =  12. 5u.  These  combine  to  yield  the  simple  result  that 

^absorption  ’  regardless  of  the  material  properties. 

The  energy  intensity  e-foiding  depth  is  given  by 


f*  8 

Sinclair's  handbook  of  aerosols  states  that 
optical  interference  between  particles  would  be  expected  if 
the  particles  were  less  than  five  diameters  apart.  Since 
the  average  spacing  between  particles  is  approximately  510 
microns  and  the  average  size  particle  is  25  microns  in  di¬ 
ameter,  the  particles  are  spaced  about  20  diameters  apart. 
Therefore,  it  is  not  necessary  to  consider  interference 

between  particles.  The  experiments  of  Churchill,  Clark,  and 
6  9 

Sliepcevich  *  enforce  this  conclusion,  since  they  deter¬ 
mined  that  the  upper  limit  of  the  ratio  of  spacing/diameter 
of  the  particles  in  a  cloud  is  more  like  1.7,  rather  than  5, 
and  any  ratio  above  1.7  indicates  that  optical  interference 
can  be  neglected. 
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The  critical  dust  cloud  density  for  consideration 
of  optical  interference  is  determined  by  using  the  criterion 


where  N  is  the  concentration  of  particles  per  unit  volume 
and  (  is  the  center-to-center  distance  between  particles. 
This  distance  is  given  by 


Since 


> 


where  P  is  the  cloud  density,  oc  is  the  dust  particle 
density,  and  r  is  the  particle  radius, 


o  -  11 

p  "  72 


p„  =  0.1 
c 


3 

In  the  special  case  where  P  =  60  lb/ft  ,  the  cloud  density 

o  C 

must  be  6  lb/ft^  for  optical  interference  effects  to  become 
significant. 

6.5  IITRI  Air -Temperature  Model 

6.5.1  General 


The  air-temperature  model  considers  the  exchange 
of  thermal  energy  between  the  dust  and  the  air  as  the  dust 
is  transported  vertically  through  the  air.  The  horizontal 
and  vertical  motion  of  the  dust  relative  to  the  air  has 
been  defined  in  the  previous  chapter.  In  the  air-temperature 
model  however,  it  is  important  to  emphasize  the  relative 
aspects  of  the  vertical  motion  since  the  air  will  be  permit¬ 
ted  to  undergo  a  change  in  volume  and,  hence,  will  experience 
some  gross  vertical  motion.  Results  of  the  dust  study, 
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together  with  soil  surveys  at  potential  sites,  indicate 
(1)  that  the  bulk  of  the  soil  mass  is  composed  of  fine  par¬ 
ticles  corresponding  to  the  finer  three  or  four  size  groups, 
and  (2)  that  the  behavior  of  these  fine  particles  during 
the  positive  phase  duration  portion  of  the  blast  wave  are 
essentially  the  same  (see  Figure  5.13).  The  behavior  of 
the  particles  of  these  size  groups  is  nearly  identical 
because  the  influence  of  gravitational  forces  (as  measured 
by  the  magnitude  of  the  terminal  velocity)  is  small.  The 
air-temperature  model,  therefore,  neglects  the  size  effects 
of  the  particles  as  far  as  the  transport  aspects  of  the  dust 
cloud  are  concerned. 

The  exchange  of  energy  between  a  particle  of  soil 
and  its  surroundings,  including  the  air  in  which  it  is  placed 
will  be  accomplished  by  the  mechanisms  of  convection  and 
radiation.  The  transfer  of  the  stored  thermal  energy  in 
the  soil  particle  will  be  governed  by  the  mechanism  of  conduc 
tion.  An  examination  of  the  conduction  problems  indicates 
that  for  the  very  small  particle  sizes  of  interest,  the  major 
ity  of  stored  energy  can  be  extracted  from  the  interior  of 
the  particle  in  a  time  of  the  order  of  one  microsecond. 

Thus,  the  conduction  process  is  not  a  limiting  process.  An 
examination  of  a  simple  radiation  problem  indicates  that  the 
loss  of  a  significant  amount  of  the  stored  energy  by  this 
mechanism  *will  require  approximately  10  sec.  The  dust  will 
be  transported  a  significant  distance  through  the  air  in 
this  period  of  time;  hence,  a  relatively  small  temperature 
change  will  occur  in  the  dust  due  to  thermal  radiation. 

The  dust  particles  are  carried  through  the  air  by 
the  local  turbulent  action  of  the  air;  hence,  very  little 
relative  motion  will  exist  between  the  air  and  the  dust 
particles.  The  convection  process  will,  therefore,  be  that 
of  natural  or  free  convection.  Normally,  radiative  cooling 
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represents  a  significant  portion  of  the  energy  transport  in 
natural  convection  heat  transfer  problems,  except  for  very 
small  objects  as  is  the  case  herein.  Estimates  of  the  heat 
transfer  coefficients  for  the  current  problem  indicate  that 
the  majority  of  the  thermal  energy  stored  in  the  soil 
particles  will  be  transferred  to  the  surrounding  air  in  times 
of  the  order  of  several  tens  of  milliseconds.  Thus,  the 
controlling  heat  transfer  process  is  that  of  natural  convec¬ 
tion.  The  air-temperature  model  is  based  upon  this  heat 
transfer  mechanism. 

The  air-temperature  model  consists  of  a  simple 
periodic  thermal  energy  balance  for  a  local  mass  of  air.  The 
dust  which  mixes  with  this  mass  of  air  approaches  thermal 
equilibrium  with  the  air  in  some  time  interval,  the  time 
interval  being  related  to  the  heat  transfer  coefficient. 

The  dust  particles  will  generally  enter  the  region  or  cell 
occupied  by  the  given  mass  of  air  at  a  temperature  which  dif¬ 
fers  from  the  current  air  temperature.  The  dust  temperature 
will  be  determined  from  the  air  temperature  of  the  cell  from 
which  it  just  left  or,  if  the  air  is  adjacent  to  the  soil 
surface,  from  the  temperature  of  the  soil  which  is  currently 
being  eroded  by  the  blast  wind.  The  soil  particles  in  a 
particular  cell  will  also  absorb  some  incident  thermal  radi¬ 
ation  which  is  continually  being  emitted  from  the  fireball. 
Since  the  vertical  velocity  of  the  dust  particles  is  small 
when  compared  to  the  acoustic  velocity  of  the  air,  any  pres¬ 
sure  waves  which  could  be  generated  would  have  propagated 
considerable  distances  compared  to  the  cell  size.  Thus,  the 
energy  balance  occurs  at  constant  pressure,  ra  her  than  at 
constant  volume.  The  transient  pressure  field  which  exists 
in  the  blast  wave  is  assumed  to  be  unaltered  by  the  existence 
of  the  dust  and  thermal  energy.  The  mass  of  air  or  the  cell 
will  undergo  a  change  in  volume  due  to  the  changing  pressure 
field  and  the  change  in  the  thermal  content  of  the  cell. 


The  cell  will  therefore  perform  work  on  its  surroundings. 

The  initial  conditions  of  the  air  at  the  time  of  arrival  of 
the  shock  wave  will  correspond  to  the  shock  state  at  the  range 
of  interest. 

6.5.2  Description  of  the  Air -Temperature  Model 

A  packet  of  hot  dust,  which  is  eroded  at  a  particu¬ 
lar  time,  enters  a  given  column  of  air.  The  horizontal  mo¬ 
tion  of  those  specific  dust  particles  is  identical  to  that 
of  the  air  column.  The  vertical  motion  is  that  which  is 
described  in  the  cold  dust  model.  The  vertical  temperature 
distribution  in  the  column  of  air  will  vary  with  time  as 
the  column  sweeps  over  the  ground  surface.  The  vertical 
temperature  distribution  at  the  site  is  obtained  by  consid¬ 
ering  the  vertical  temperature  distribution  in  the  column 
of  air  which  is  currently  at  the  site  location.  The  model, 
therefore,  considers  the  transient  temperature  variations  of 
a  number  of  air  columns  up  until  the  time  that  they  sweep 
past  the  site.  Time  measured  along  the  path  of  the  air  column 
is  related  to  the  observation  at  the  site  by  the  equation 


(6.22) 


where 

t  =  time  following  air  column, 

t  =0  when  air  column  is  set  into  motion, 
tQbs  “  binie  at  site  location, 
tobs  =  0  when  shock  arrives  at  site, 

R  =  current  displacement  of  air  column,  and 
U  =  shock  velocity. 

For  a  uniform  ground  condition,  all  of  the  air  columns  experi 
ence  the  same  temperature  transient  except  for  a  time  shift. 
The  solution  to  this  simple  case  is  then  obtained  by  consid¬ 
ering  the  temperature  variation  for  a  single  air  column  as 
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it  moves  across  the  ground,  and  then  connecting  the  time 
scale  in  accordance  with  Equation  6.22. 

The  column  of  air  is  divided  up  into  a  number  of 
Cells  whose  initial  height  is  6y.  The  weight  of  air  per 
unit  ground  area  is 

w  =  pga6y  (6.23) 

where 

w  =  weight  of  air,  and 

=  air  density  (shock  density). 

The  energy  balance  for  each  cell  is  made  periodically  at  an 
interval  of  time  needed  to  approach  thermal  equilibrium 
between  the  dust  and  the  air  in  the  cell.  This  time, 
t  ,  is  determined  from  the  heat  transfer  coefficient  and 
represents  a  variable  for  the  current  problem. 

The  maximum  vertical  velocity  of  the  dust  is  KVq 
where  K  is  the  transport  coefficient  and  VQ  is  the  peak  air 
velocity.  It  will  be  convenient  to  select  an  initial  cell 
size  DY  such  that  the  dust  particles  will  not  move  a  distance 
greater  than  6y  during  a  t  .me  step  t  .  Furthermore,  it  will 
be  desirable  from  a  computational  point  of  view  to  require 
that  all  of  the  cells  have  the  same  initial  size  and,  there¬ 
fore,  the  same  weight  of  air.  Thus, 

6y  =  te  K  V0  .  (6.24) 

Recalling  that ,  in  the  absence  of  gravitatioi  al 
effects,  the  cloud  height  H  is  governed  by  the  differential 
equation 


(6.25) 
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where 

t  =  time,  and 

V  =  horizontal  wind  velocity, 


the  current  model  will  be  limited  to  the  positive  phase  dura¬ 
tion  of  the  blast -induced  winds;  hence,  the  absolute  velocity 
restriction  can  be  dropped  from  Equation  6.25.  It  will  be 
convenient  to  select  a  somewhat  nonuniform  time  step  At  in 
order  that  the  change  of  the  height  of  the  cloud  during 
each  time  step  be  uniform.  During  the  early  portion  of  the 
dust  cloud  growth,  when  the  horizontal  wind  velocity  is 
changing  slowly,  the  time  step  will  be  nearly  uniform. 

The  initial  time  step  is  equal  to  tg;  hence,  the  change  in 
height  of  the  cloud  during  each  time  step  is  equal  to  the 
initial  cell  height  6y.  It  should  be  noted  that  these 
conditions  apply  in  a  Lagrangian  sense,  i.e.,  in  the  absence 
of  any  vertical  dilatation  of  the  air  column.  The  actual 
height  of  the  dust  cloud  would  depend  upon  the  temperature 
distribution  within  the  column  of  air. 

A  Lagrangian  cell  configuration  can  now  be  estab¬ 
lished.  This  configuration  is  illustrated  in  Figure  6.16(a). 
The  independent  variables  are  the  Lagrangian  height  y^  and 
the  time  t^,  where  the  indicies  i  and  j  refer,  respectively, 
to  a  specific  time  and  a  specific  cell.  The  cell  is  defined 
by  the  coordinate  of  its  upper  edge.  The  top  of  the  cloud 
is  thus  given  by 


Hi  =  (i-1)  6  y 


and  the  time  is  given  by 


(6.26) 


(6.27) 
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|  is  the  average  value  of  the  horizontal  air  velocity  during 

Ithe  time  step. 

The  horizontal  displacement  of  the  air  column 
!  during  a  time  step  is 

DR  =  V.  (t^  -  t^ 

and  is  therefore  a  constant,  equal  to 

DR  =  .  (6.28) 

Furthermore,  the  quantity  of  dust  picked  up  per  unit  area 
during  a  time  step  is 

W  =  Fap gVi  <t.  -  t-.p 

where 

a  *  erosion  coefficient, 

P  *  bulk  density  of  soil, 

b 

F  =  land  use  factor, 
and  is  therefore  a  constant,  equal  to 

P  6y 

W  =  Fa  .  (6.29) 

The  total  horizontal  displacement  of  the  air  column  is 


K±  =  (i-1)  DR  =  Hj/K; 


(6.30) 


li.nce,  the  observation  time  ^  defined  in  Equation  6.22 

can  be  rewritten  as 


t-  =  t  -  (i-PDR 
tobs,i  ti  U 


(6.31) 
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A  constant -pressure  energy  balance  is  made  for 
each  cell  at  each  interval  of  time  in  accordance  with  the 
illustration  given  in  Figure  6.16(b),  wn>.re  is  the  quan 
tity  of  dust  per  unit  area  in  the  cell  at  the  iEll  time. 
Similarly,  w!  .  and  wV.  are  the  quantities  of  dust  per  unit 
area  which,  respectively,  leave  the  j-1  cell  and  enter  the 
j—  cell,  and  which  leave  the  j~  cell  and  enter  the  j+1 
cell  during  the  i—  time  interval.  Clearly, 


,TU  _ 

W. .  —  W. 


i,j+l  * 


Therefore,  we  need  only  deal  with  the  quantities  w^^  and 

wl^.  A  simple  formula  for  w^  can  now  be  developed  by 

considering  the  contribution  of  dust  in  the  cell  from  all  of 

ttl 

the  dust  packets.  The  total  quantity  of  dust  in  the  k— 

packet  is  w  and  is  distributed  to  its  current  height  H  in  a 

linear  manner.  This  is  illustrated  in  Figure  6.17.  The 

th 

quantity  of  dust  in  the  j —  cell  is 


v;ijk  =  Trl1  -  iW  • 


But , 


and 


H  =  6  y  (i-k), 
y  =  6y  (j  -  |) 


dy  =  6  y  ; 


hence , 


2W  ,  ,  3 

iik - 7  i-K-j  +  7  . 

1JR  (i-k)Z  l  J 


(6.32) 


For  a  uniform  soil  condition  the  total  quantity  of  dust  in 
the  cell  w^j  is  given  by  summing  over  all  of  the  k—  dust 
packets.  (It  should  be  noted  that  dust  from  packets  which 
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have  entered  the  overall  cloud,  but  which  cannot  reach  the 

th  iVi 

j —  ceil  at  the  i~—  time,  must  be  excluded.)  Hence, 


i-j+1 

r- 

w.  .  *  >  w. 

k=l 


(6.33) 


or 


wij  "  wi-l,j  +  wijl  ’ 


which  becomes 


wtj  -  wt.t  j  +  Ml CL+JLQ.  . 
1J  1  i,J  (i-l)Z 


(6.34) 


th 


The  quantify  of  dust  per  unit  area  from  the  k — 

dust  packet  which  enters  the  j —  cell  during  the  i=—  time 

interval  is  w^^.  Figure  6.18  illustrates  the  distribution 

of  the  dust  from  this  packet  as  well  as  the  bounds  of  the 

t~h 

subgroup  of  particles  which  enters  the  j —  cell.  The  dust 

moves  in  such  a  manner  that  the  quantity  of  dust,  bounded  by 

th 

any  pair  of  rays  passing  through  the  k—  origin,  is  constant. 

Vi  Vi 

Therefore,  the  dust  entering  the  j —  cell  during  the  i— 
time  is  bounded  by  the  ray  passing  through  the  point  (i-l,j-l) 
and  by  the  ray  passing  through  the  point  (i,j-l).  The  fol¬ 
lowing  equations  result: 


wijk 


2W  j 

l 


M 


dy 


but 

and 


H  a  6y  (i-k)  *  H. 


i  _ 


y'“  yj-i  Bik 


ik 


where 


Bik  ”  h 


Hik 


i-l,k 
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and 


y  =  \  (yi-i  +  y')  . 

dy  =  y'  -  ; 


hence 


(-IWl.  «■?» 


For  a  uniform  soil  condition  the  total  quantity  of  dust 

a.L  i.t. 

entering  the  j—  cell  during  the  i— -  time  interval  is  given 
by  summing  all  of  the  k— -  oust  packets,  noting  again  that 
some  packets  must  be  excluded.  This  quantity,  then,  is  given 
by 

i-j 

L  wijk  <6-36> 

k=l 


w!  .  = 


ij 


or 


w.'  .  -  w!  ,  .  +  w! . 

ij  i-l. J  ijl 


which  becomes 


(6.37) 

It  should  be  further  noted  that  both  w.  .  and  w.  ..  are  zero  in 

J-J  l-jK 

the  region  above  the  dust  cloud  (j-i<0). 


Two  types  of  nonuniform  or  cutoff  ground  conditions 
have  also  been  treated.  These  are  the  situations  where 
(1)  erosion  stops  because  of  an  erosion  cutoff  depth  z  has 
been  reached,  or  (2)  the  air  column  has  reached  a  range  where 
the  land  use  factor  has  changed.  These  two  problems  are 
similar  if  we  determine  two  different  values  for  the  quantity 
of  soil,  W,  in  a  dust  packet.  Let  these  be  and  V^,  where 
*  0  for  the  case  of  no  erosion.  The  corresponding  cell 
configuration  is  illustrated  in  Figure  6.19.  Regions  0  and 
2  are  those  treated  in  the  preceding  paragraph.  Region  1  is 
the  new  region  introduced  by  the  change  in  the  land  use 


p 
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factor  or  by  erosion  depth  cutoff.  The  boundary  between 
Regions  1  and  2  is  defined  by  the  index  i  «  ic.  If  erosion 
stops  after  a  depth  z  has  been  reached,  then 

ic  =  ^  (6.38) 


or  the  land  use  factor  changes  at  a  range  R  in  front  of  the 
site.  Then, 


(6.39) 


since  at  the  i-—  time,  the  column  is  (by  definition)  at  the 
site. 

For  cells  in  Region  1,  the  quantity  of  dust  w^ 
in  the  cell  is  determined  as  before;  however,  the  summation 
is  divided  into  two  group;?:  one  involving  W^,  and  the  other 
involving  W2.  The  results  are 


+ 


where  w^  is  the  contribution  from  zone  1  (W  =  W^)  and  w^ 
is  the  contribution  from  zone  2  (W  =  W2) .  Therefore, 


Wi-l,j 


+ 


wijl 


w. 


1J,1C 


w.  .  . 

1J,1C 


(6.40) 


or 


2W, 


w. 


.  +  - 1— 

1'1,J  (i-l)* 


(i-j+1/2) 


2(W2-W1) 

,2 

(l-ic) 


(i-ic-j+1/2) , 

(6.41) 


recalling  that  for  the  erosion  cutoff  case,  *  0. 

A  similar  treatment  for  w^  results  in 


w! .  =  w!  .  .  +  w!  ..  -  w! .  .  +  w! .  . 

1J  Ijl  1J,1C  1J,LC  » 


(6.42) 
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or 


2W2(j-2)  (*  n-2)  a-3/2)l 

wij  "  Wi-I,j  +  (i-l)(T-TF  ) 1  ■  ( 

2(v2-v{; 


(i-ic) (i-ic-1) 


H-2)(i-ic-l/2)l 
-  *  (i-ic)  (i-ic-Tjy  * 


(6.43) 


The  energy  balance  on  the  cell,  as  illustrated  in 
Figure  6.16,  is 

Eij  +  Eij  ‘  Eij  ■  Eij  =  +  cwij)(Tij  ' 


The  left-hand  side  of  the  equation  represents  the  net  energy 
flux  which  consists  of  the  energy  convected  into  the  cell  by 
the  dust  w! .  and  is 


Wij 


c  T. 


•  •  t 


where  c  is  the  specific  heat  of  the  dust,  and  T^^_^  is  the 
air  temperature  in  the  j-1  cell  at  t±.  The  energy  deposited 
in  the  cell  by  direct  thermal  radiation  is  ,  and  is  ob¬ 
tained  from 


Eij  "  (Qij  ‘  Vj-l)(ti  '  ti-l) 


(6.45) 


where 


=  Qn  exP<"  wij/wo  cos 


(6.46) 


and 

0  =  thermal  flux  incident  upon  the  cloud, 

o 

W*.  =  quantity  of  dust  above  the  j —  cell  at 
1-j  time  =  t^, 

w  =  weight  absorption  coefficient  permit  area, 
°  and 

9  =  angle  of  incidence  of  radiation. 
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The  energy  convected  rut  of  the  cell  by  the  dust  leaving  the 
cell  is 


w'.' . 


c  T  .  ,  .  , , 


or 


Ei,j+1  Wi,j+1  c  Ti-l,j+l 


(6.47) 


The  expansion  work  that  the  cell  does  upon  its  surroundings 
is  eV.  and  is  determined  from 


P 


i-1 


) 


(6.48) 


where 

=  pressure  at  t  =  t.  ,  and 

t-  h 

S. .=  height  of  j —  cell  at  t  -  t.  . 
ij  1 


The  deviation  of  S.j  will  be  presented  in  the  following 
paragraph.  The  right-hand  side  of  Equation  6.44  is  the 
change  in  the  internal  energy  of  the  cell.  The  mass  of  air 
in  the  cell  is  w,  whereas  c  is  the  specific  heat  of  air  at 
constant  pressure.  T. -  -  T.  1  .is,  of  course,  the  tempera 
ture  change. 

The  air  cells  will  undergo  a  change  in  volume  due 
to  both  the  thermal  input  and  the  changing  pressure  field. 
In  the  absence  of  dust,  the  volume  change  will  only  be  due 
to  the  changing  pressure  field.  It  is  assumed  that  this 
volume  change  will  occur  in  the  horizontal  direction.  The 
remainder  of  the  volume  change  will  occur  in  the  vertical 
direction  and  v/ill  transport  the  air  vertically  (see  Fig¬ 
ure  6.20).  The  volume  v  of  a  cell,  per  unit  width,  is 


v  =  SL 
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Figure  6.20  Cell  Dilation 


where 

S  =  cell  height  and 
L  =  cell  length. 

Differentiating,  we  obtain 

dv  _  dS  .  dL 

v  'S-  +  IT  • 

For  an  isentropic  expansion  (the  cause  of  the  horizontal 
dilitation) , 


dL  _  v  dP 

r- -c;r 

where  c  and  c  are  the  specific  heats  at  constant  volume 
v  p 

and  at  constant  pressure,  respectively,  and  where  P  is  the 
pressure.  Thus,  the  vertical  dilitation  must  be 


dS  dv  dL 
3  v  XT 


But  for  a  perfect  gas. 


dv  dT  dP 
v~  ’  T  '  TT  • 


Hence , 


dS  dT  ,  v  1 1  dP 

r  r  \  v  r  ' 


(6.49) 


Equation  6.49  can  now  be  rewritten  as 


c  -  c  I  i  dT  L  Sr  dP 

si  '  si-l  1  +  T"  *  7~  • 


(6.50) 


The  center  height  of  the  cell  (Y^^j)  is  given  by  surnning 
over  the  appropriate  cells  as 
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(6.51) 


+  7 


+  S 


The  energy  balance  Is  then  performed  sequentially 
on  each  cell  time  or  ij  combination.  The  number  of  cells 
which  must  be  treated  depends  upon  the  cell  size  and  the 
overheight  for  which  information  is  required.  Two  specific 
programs  were  written  in  FORTRAN  IV  based  upon  the  above 
model. 

One  program  deals  with  a  uniform  ground  condition 
and  includes  erosion  cutoff.  This  program  follows  a  single 
air  column  and  then  adjusts  the  time  scale  to  the  observa¬ 
tion  time  at  the  site.  The  temperature  is  computed  at  some 
fixed  elevation  Y  by  assuming  a  linear  temperature  variation 
between  adjacent  cell  centers.  The  dust  densities  are 
for  each  cell  and  also  are  computed,  since  both  the  quantity 
of  the  dust  in  each  cell  and  the  volume  of  each  cell  is  known. 
For  example. 


(6.52) 


The  dust  density  at  some  height  Y  is  obtained  by  assuming 
a  linear  variation  between  cell  centers. 

The  second  program  treats  a  biuniform  ground  condi¬ 
tion,  i.e.,  one  in  which  two  uniform  ground  conditions  exist. 
The  boundary  between  these  zones  is  located  a  distance  RC 
in  front  of  the  site.  This  code  requires  a  considerably 
greater  computational  time  since  the  temperature  variation 
of  many  air  columns  must  be  evaluated.  Some  savings  in 
computational  time  can  be  achieved  by  noting  that  (1)  during 
the  early  portion  of  the  problem  the  air  columns  passing 
the  site  have  only  swept  past  the  one  zone,  hence,  the  uni¬ 
form  ground  approach  can  be  used,  and  (2)  that  a  first 
portion  of  each  air  column  solution  in  the  two  zone  case  is 
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identical  to  the  preceding  air  column  problem.  These  pro¬ 
grams  are  presented  in  Appendix  A.  Some  typical  results 
which  illustrate  the  nature  of  this  phenomenon  are  presented 
in  Figures  6.21  and  6.22.  These  figures  show  the  tempera¬ 
ture  distribution  in  the  air  at  various  times.  Initially, 
all  of  the  air  is  heated  to  a  uniform  temperature  by  the 
passage  of  the  shock  wave.  Then,  as  the  hot  soil  is  eroded 
and  enters  the  air  columi.,  the  air  temperature  over  a  portion 
of  the  column  height  increases  significantly  (see  t  =  0.02 
and  t  =  0.13  sec).  However,  as  the  process  continues,  all 
of  the  hot  soil  has  been  eroded  and  the  cooler  dust  begins 
to  lower  the  temperature  of  the  bottom  portion  of  the  air 
column  (see  t  =  0.27  and  0.38  sec).  Finally,  at  later  times, 
the  unheated  dust  (at  the  ambient  temperature)  enters  the 
dust  cloud  and  further  cools  the  air.  A  very  definite 
cooling  or  quenching  phase  does  exist.  The  maximum  air 
temperature  at  any  height  decreases  with  increasing  height 
as  indicated  in  Figure  6.21. 

The  results  illustrated  in  Figure  6.22  are  for  the 
case  where  cutoff  occurs  at  a  time  t  *  0.39  sec.  Once  cutoff 
occurs  the  temperature  at  the  lower  portion  of  the  air  column 
begins  to  stabilize  and  the  upward  motion  of  the  thermal 
pulse  is  considerably  slowed.  A  comparison  is  made  with  the 
noncutoff  case  at  t  =  135  sec  and  indicates  that  a  consider¬ 
able  difference  occurs.  The  temperature  distribution  at 
t  =  4.24  sec  indicates  that  the  heated  air  is  undergoing 
a  small  change.  It  is  important  to  note  that  the  air- 
temperature  model  used  herein  assumes  that  the  heated  layer 
of  air  will  not  become  unstable  and  break  up  during  the  time 
of  interest.  Obviously,  the  time  required  for  the  layer  to 
become  unstable  will  depend  upon  both  temperature  intensity 
and  the  magnitude  of  the  horizontal  and  vertical  temperature 
gradient's  in  the  air. 


fee 


Figure  6.1']  Vertical  Temperature  Distribution  in  Air  Column 


Height,  Y,  feet 


0  1000 
Air  Temperature,  T,  °F 


Figure  6.22  Vertical  Temperature  Distribution  in  Air 
Column  with  Erosion  Cut-off 
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The  temperature  variation  at  a  specific  elevation 
is  presented  in  Figure  6.23  for  a  typical  problem.  This 
figure  presents  the  results  for  a  variety  of  all  sizes,  or 
what  is  equivalent  to  a  variety  of  heat  transfer  coefficients, 
as  well  as  the  temperature  variation  for  the  case  of  no  dust. 
The  air  first  experiences  the  shock  heating  and  then  the 
gradual  cooling  as  the  pressure  field  decays.  When  the  dust 
reaches  the  elevation  of  interest  the  air  temperature  in¬ 
creases  rapidly,  reaching  a  peak,  and  then  decays  more  gradu¬ 
ally.  The  value  of  the  heat  transfer  coefficient  is  not  well 
defined  and  will  be  discussed  in  more  detail  in  a  later  sec¬ 
tion  of  this  chapter;  however,  it  should  be  noted  from  the 
results  presented  in  Figure  6.23  that  this  coefficient  is  an 
important  variable --one  which  has  a  considerable  amount  of 
influence  upon  the  peak  temperature.  In  particular,  the 
peak  temperature  is  not  approaching  a  stable  value  as  the 
cell  size  is  reduced. 

6.5.3  Analytic  Solution 

As  a  result  of  the  sensitivity  of  the  peak  air 
temperature  to  the  value  of  the  heat  transfer  coefficient  or 
to  the  cell  size,  an  investigation  was  made  in  an  attempt 
to  study  the  nature  of  the  transport  of  the  energy  in  the 
dust  cloud.  A  simple  problem  was  formulated  in  which  all 
of  the  secondary  effects  were  omitted.  The  problem  studied 
was  that  of  an  impulsively  induced  wind  eroding  a  uniform 
soil.  The  wind  velocity  is  constant  at  a  value  of  VQ,  the 
soil  temperature  is  uniform  with  depth  at  a  value  of  TQ,  and 
the  initial  air  temperature  is  zero.  Such  a  simple  problem 
does  not  involve  a  characteristic  dimension.  Hence,  it  would 
be  reasonable  to  expect  that  a  similarity  solution  exists  in 
which  the  air  temperature  is  a  function  of  a  similarity 
variable,  as 

T(y,t)  -  T  (y/t)  (6.53) 
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Time,  t,  seconds 

Figure  6,23  Influence  of  Cell  Size 
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where 

y  =  Lagrangian  height  above  surface  and 
t  =  time. 


Recall  that  the  dust  density  is  given  by  the  similarity  solu 


tion 


ap 


p(c)  =  -JC2-  [c-  1  -  lr'C  ] 


(6.54) 


where 

C  =  y/t  . 

Consider  the  differential  element  dy  (illustrated 
in  Figure  6.24).  The  net  mass  flux  of  dust  into  the  cell 
equals  the  rate  of  increase  of  dust  in  the  cell;  for  example 

—  dy  =  -  (6.55) 

sy  y  at 


where 

w  =  weight  of  dust  and 
w1  =  weight  flux  of  dust. 

Similarly,  an  energy  balance  yields 

-  ew'  dy  -  o  X  dy  =  |f  ,  (6.56) 

where 

E  =  (wc  +  wc)  T,  (6.57) 

and 

w  =  weight  of  air  in  differential  element,  and 
c  =  specific  heat  at  constant  pressure  of  air. 


Using  Equations  6.55  and  6.57,  we  can  rewrite  Eauation  6.56 
as 


ST 

St 


cw 


cw  +  cw 


) 
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igure  6.24  Differential  Element 


since 


w  =  pdy  and 
w  *  pdy; 

where 

p  «  dust  density  and 
o  =  air  density. 

The  energy  equation  then  becomes 


ST 

St 


cw 


Pc  +  Pc 


ST 

ay 


(6.58) 


The  quantity  w'  can  be  derived  from  the  dust  density 
distribution  illustrated  in  Figure  6.25.  The  weight  of  dust 
above  a  height  y^  is 


(6.59) 


which  upon  integration  becomes 

•2  i 


2HP  a 

W1 - JT~ 


1  ^  1 

1  *  T  +  C1  ln  C1 


The  total  amount  of  dust  raised  WQ  is 

HP  a 

=  — p-—  =  V  to  a  . 
o  K  OS 


(6.60) 


(6.61) 


However,  since 


dW, 

w' =  -at-  • 


we  obtain 


w*  =  VQDsa  (1  - C  ) 


(6.62) 


where  has  been  set  equal  to  Q  .  The  differential  form  of 
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the  energy  equation  (Equation  6.58)  can  now  be  rewritten  as 


f£  -  -  Wo G<C>  § 

(6.63) 

where 

2 

(6.64a) 

*  B+2(C  '  1  "  ln  C  ) 

and 

(6.64b) 

r  -  2SJL 

B  P„ca 

Equation  6.63  does  not  satisfy  the  conditions  for  a  similarity 
solution  except  ir.  the  trivial  case  where  both  *T /at  and 
*T/dy  vanish.  Therefore  we  must  consider  the  possibility  of 
a  discontinuous  temperature  field.  This  idealized  problem 
was  programmed;  however,  since  the  computation  starts  with 
a  single  cell,  the  temperature  distribution  is  initially 
uniform.  We  can  only  expect  that  the  similarity  solution 
would  evolve  from  this  starting  procedure.  The  results  of 
this  calculation  are  presented  in  Figure  6.26  and  clearly 
demonstrate  that,  as  the  number  of  time  steps  N  becomes 
large  the  temperature  distribution  approaches  a  dir  continu¬ 
ity  at  the  normalized  height  ^  *  0.9  for  this  particular 

example . 

If  we  examine  the  function  G(C)  we  note  that  it  has 
a  maximum  at  approximately  C  *  C i  (see  Figure  6.28),  The 
function  CO  is  the  ratio  of  energy  influx  to  the  cell  to 
the  thermal  capacity  of  the  cell.  Differentiating  Equa¬ 
tion  6.64a  results  in  a  relation  for  the  location  of  the 
temperature  discontinuity  £  *  Cc 

1  »  BCc  +  ^c  '  2Cc  ln  Cc*  (6,65) 

The  similarity  solution  is. illustrated  in  Figure  6.27. 
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The  above  solution  can  be  verified  by  applying 
the  overall  energy  equation  to  the  system.  That  is, 


:  wA  »  / 

00  il 


cP  (C)  +  cP  T(C)  dc  . 


By  letting 


T(C)  *  T 


0  >  C  >  C 


(6.66) 


and 


T(C )  *  0 


Cc  -  C  —  *> 


Equation  6.65  is  confirmed. 

The  air  temperature  can  be  expected  to  reach  or 
approach  the  maximum  soil  temperature  if  the  heat  transfer 
coefficient  is  very  high.  Thus,  the  results  illustrated 
in  Figure  C . 23  are  confirmed  and  the  limits  have  been 

established. 

6.5.4  Heat  Transfer  Coefficient 

The  air  temperature  model  contains  only  one  coef- 
ficient  which  must  be  evaluated.  This  is  the  heat  transfer 
coefficient  for  natural  convection  between  the  small  irregu¬ 
larly  shaped  soil  particles  and  the  surrounding  air .  This 
coefficient  is  used  to  evaluate  an  energy  exchange  response 
time  t  which,  In  turn,  is  used  to  evaluate  a  cell  size. 

There  are  other  coefficients  involved  in  the  overall  model . 
These  coefficients  h.  ve  all  been  discussed  and  evaluated  in 

the  preceding  chapter. 

Heat  transfer  data  which  are  pertinent  to  thin 
problem  are  both  sparse  and  unreliable.  Data  exist  on 
simple  shapes  such  as  spheres  and  horizontal  cylinders;  • 
however,  these  data  are  generally  restricted  to  oojteis  -.hose 
dimensions  are  of  the  order  of  inches  in  size.  There  are  some 

data  on  fine  horizontal  wires. 
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The  natural  convection  phenomenon  can  be  analyzed 
in  terms  of  three  dimensi  nless  groups.  These  are: 
Nusselt's  Numbf  *  hD/k, 


and 

where 


Grashof's  Numuer  *  D^o^ST/u^, 
Prandtl's  Number  *  cu/k 


D  *  diameter, 

h  *  heat  transfer  coefficient, 
k  *  themal  conductivity  of  air, 
p  *  density  of  air, 
u  *  viscosity  of  air, 

3  *  coefficient  of  thermal  expension  for  air, 
T  *  temperature  difference,  and 
c  *  specific  heat  of  air 


Experimental  data  for  horizontal  cylinders  (and  wires)  have 
been  obtained  for  a  variety  of  fluids  and  fluid  properties. 
Thus,  the  interrelationship  between  these  three  dimensionless 
groups  has  been  established  for  most  engineering  applications. 


For  objects  of  the  order  of  inches  in  size,  the 
heat  transfer  coefficient  is  proportional  to  the  one-fourth 
root  of  the  size.  As  the  size  decreases,  however,  the  de¬ 
pendence  approaches  an  inverse  proportionality.  Due  to  the 
strong  dependence  of  the  Grashof  number  on  the  diameter,  the 
current  problem  falls  in  a  range  where  little  data  exist. 
Because  of  this  lack  of  pertinent  data,  an  experimental  pro¬ 
gram  was  designed  and  conducted  to  determine  the  dependence 
cJ:  the  heat  transfer  coefficient  for  natural  convection  on 
the  size  of  typical  soil  particles.  The  results  of  this 
experiment  are  reported  in  Appendix  E. 


The  value  of  the  heat  transfer  coefficient  as  a 
function  of  particle  size  is  presented  in  Figure  6.29.  The 
smaller  the  particle  size,  the  larger  the  value  of  the  heat 
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transfer  coefficient.  If  a  mixture  of  particle  sixes  was 
used  then  the  effective  value  of  the  heat  transfer  coef¬ 
ficient  appeared  to  be  chat  associated  with  the  average 
particle  size.  Since,  in  general,  soils  are  composed  of 
a  mixture  of  particle  sizes,  it  was  necessary  to  examine 
the  size  distribution  curve  and  select  some  mean  size. 

For  many  of  the  typical  soils  of  interest  this  mean  size 
is  approximately  20  microns  in  diameter.  The  value  of  the 
heat  transfer  coefficient  for  this  size  is  420  Btu/ft^- 
hr-°F . 


The  energy  exchange  or  heat  transfer  response 
time  is  related  to  the  heat  transfer  coefficient  by  the 
relation: 


0  c  D 
=  •>  c  p 

3C5HTI 


(6.67) 


The  term  enclosed  in  the  brackets  is  the  exponential  time 
constant  associated  with  the  transient  heat  transfer  pro¬ 
cess  and  the  factor  3  is  a  multiplier  that  corresponds  to 
the  transfer  of  95  percent  of  the  thermal  energy  from  the 
soil  particle  to  the  air.  The  95  percent  criterion  is  an 
arbitrary  criterion  that  could  be  modified  if  desired.  The 
heat  transfer  response  time  corresponding  to  the  20-micron 
size  particle  was  approximate  ljT  8.5  milliseconds.  A  nominal 
value  of  10  milliseconds  is  suggested  for  the  heat  transfer 
response  time  for  those  problems  where  special  soil  condi¬ 
tions  do  not  exist. 


6.5.5  Stability  of  the  Heated  Air  Layer 

The  air  near  the  ground  surface  is  heated  and  then 
cooled  oy  the  passage  of  dust  through  this  region.  In  gen¬ 
eral,  the  horizontal  temperature  gradients  are  small  and  the 
transient  temperature  field  is  one-dimensional  in  structure. 
The  heated  air  near  the  surface  is  bounded  above  and  below 
by  cooler  air  and  represents  an  unstable  or  metastable  con“ 
figuration.  The  hot  air  tends  to  rise  due  to  local  buoyancy 


* 
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forces,  although  the  migration  of  dust  through  this  region 
continuously  alters  the  temperature  distribution,  hence  the 
buoyancy  force  distribution.  The  buoyancy  forces  are  large 
initially,  but  decrease  rapidly.  It  is  necessary,  therefore, 
to  determine  whether  the  predicted  temperature  fields  are 
stable  or  whether  the  buoyancy  forces  overturn  the  unstable 
air  locally  and  grossly  change  the  temperature  distribution. 
This  subsection  discusses  the  buoyancy  problem  and  develops 
and  approximate  criterion  for  evaluating  the  stability  of 
the  heated  air  layer.  Suggestions  are  then  made  for  modi¬ 
fying  the  temperature  distribution  when  overturning  occurs. 

The  local  thermal  and  displacement  field  is  very 
complex  since  it  involves  not  only  the  exchange  of  thermal 
energy  between  the  dust  and  the  air  but  also  the  motion  of 
the  dust  or  soil  particles.  The  motion  of  the  dust  is  pos¬ 
tulated  to  move  vertically  with  respect  to  the  air  and  the 
air  is  restricted  to  move  horizontally.  The  overturning  of 
the  resultant  heated  air  layer  necessitates  the  introduction 
of  simultaneous  horizontal  and  vertical  motion  of  both  the 
dust  and  the  air,  and  greatly  increases  the  complexity  of 
any  analysis  that  can  be  used  to  evaluate  the  stability 

prob lem. 


Because  similar  stability  problems  occur  in  the 
field  of  meter ology ,  the  pertinent  literature  was  re- 
viewed.  The  convection  process  involving  the  two-dimensiona 
nonsteady  motion  of  heated  masses  of  air  and  the  fo‘Ta^Q0g  ^ 
deve lonment  of  water  vapor  clouds  have  been  reviewed  '  • 

and  appear  to  occur  at  a  variety  of  size  scales  and  involve 
a  number  of  discrete  phenomena.  Numerical  methods'- have 
been  developed  for  computing  the  motion  of  the  air  in  some 
cases.  These  computational  efforts  are  very  large  and, 
therefore,  expensive,  and  are  greatly  limited  by  frequent 
numerical  instability  problems.  These  specific  me-,  ods 
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are  not  directly  applicable  to  the  dust-ladened  heated  air 
problem  and  a  very  sizable  development  effort  would  be  nec¬ 
essary  to  establish  a  realistic  computer  program  for  the 
dust  problem. 

The  heated  air  layers  are  unquestionably  unstable 
in  many  situations  and  a  number  of  local  effects  will  always 
exist  that  will  result  in  small  pertubations  which,  in  turn, 
will  trigger  the  overturning  process.  It  is  necessary  that 
the  time  or  duration  needed  for  this  overturning  process  be 
established  and  then  compared  with  that  time  interval  during 
which  a  specific  unstable  temperature  distribution  exists. 

If  the  overturning  duration  is  greater  by  some  margin,  then 
the  preheated  air  temperature  distribution  will  change  to  a 
new  distribution  in  accordance  with  the  air  temperature 
model.  On  the  other  hand,  if  the  overturning  duration  is 
smaller,  then  we  can  expect  that  the  heated  air  layer  will 
overturn  and  mix  with  the  surrounding  air.  The  resulting 
air  temperature  distribution  will  be  much  more  uniform  and 
a  more  stable  configuration  will  result.  The  scale  of  mix¬ 
ing  cannot  be  readily  established,  however  engineering 
estimates  that  are  suitable  for  this  program  can  be  made. 

Typically,  the  vertical  air  temperature  distribu¬ 
tion  at  specific  times  will  be  similar  to  those  presented 
in  Figure  6.21  and  the  temperature  variation  at  a  nominal 
elevation  will  be  similar  to  that  presented  in  Figure  6.23. 
If  erosion  cutoff  occurs,  then  the  temperature  distribution 
will  be  much  more  uniform  over  the  lower  portion  of  the  air 
region  (Figure  6.22).  For  near-overhead  burst  situations, 
the  air  temperature  increases  much  more  slowly  than  indi¬ 
cated,  and  high  temperatures  may  persist  for  periods  be¬ 
tween  five  and  ten  seconds.  We  might  expect  that  overturn¬ 
ing  would  occur  in  these  cases. 
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where 


The  buoyancy  force  Fb  for  a  mass  of  heated  air  is 


Fb  -  p3  7 


(6.68) 


o  =  mean  air  density 

3  =  coefficient  of  expansion  *  ^  for  Perfect  Sas 

T  =  absolute  temperature 
9  =  temperature  difference. 

The  mass  involved  will,  include  both  the  mass  of  the  air  and 
the  mass  of  the  airborne  dust.  The  critical  acceleration  a 
of  this  region  of  heated  air  will  be 


_  gFb  o3g9  (6.69) 

a  p  +  od  2(p  +  P<j) 

The  temperature  difference  0  will  exist  over  some  vertical 
height  L  and  the  overturning  process  will  be  well  developed 
when  the  air  mass  has  been  displaced  this  distance.  The 
buoyancy  force  will  decay  during  the  overturning  process  and 
the  mean  acceleration  will  be  somewhat  lower  than  the  initial 
acceleration.  A  factor  of  one-half  has  been  introduced  as  a 
multiplier  of  the  initial  acceleration  to  compensate  approxi¬ 
mated  for  the  reduction  of  the  driving  force  as  well  as  for 
the  increase  in  mass  which  will  occur  due  to  air  entrainment 
and  wake  effects.  The  time  ts  required  for  the  mass  of  air 
to  be  displaced  distance  L  is  then 

•  8L(p  +  Qd)  (6.70) 

ts  ~  opge 

The  dust  and  air  density  variation  as  well  as  the  air  temper¬ 
ature  distribution  with  vertical  height  are  illustrated  in 
Figure  6.30  for  a  rather  typical  situation.  The  temperature 
difference  is  approximately  300°F  and  occurs  over  a  distance 
of  approximately  20  ft.  The  dust  and  air  density  vary  some¬ 
what  over  this  region  but  the  use  of  average  values  in  Equa¬ 
tion  (6.70)  yields  an  estimate  of  the  overturning  time  as 
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Figure  6.30  Buoyancy  Forces  in  Heated  Air  Layer 


approximately  5  seconds.  The  temperature  transient  for  this 
particular  case  lasts  for  approximately  one  second  and  is, 
therefore,  stable. 

In  those  cases  where  overturning  occurs  the  heated 
air  has  mixed  with  cooler  air,  and  the  region  originally 
filled  by  the  heated  air  is  occupied  by  air  from  the  sur¬ 
rounding  region.  The  temperature  distribution  cannot  be 
well  established;  however,  by  conserving  the  overall  energy 
of  the  system,  one  can  assume  that  the  air  temperature  in 
that  particular  region  will  be  equal  to  the  average  tempera¬ 
ture  over  a  region  two  to  three  times  as  large  as  the  orig¬ 
inal  region.  The  following  approximate  procedure  is 
suggested: 

(1)  Examine  temperature  distribution  at  time  t 
and  evaluate  overturning  time,  ts. 

(2)  Examine  temperature  distribution  at  a  time 
ts  later 

if  (a)  the  temperature  distribution  has 

changed  to  a  more  stable  distribution, 
the  process  is  stable 
if  (b)  the  temperature  distribution  is 

approximately  the  same  or  more  un¬ 
stable,  then  overturning  will  occur. 

If  so  proceed. 

(3)  Determine  average  temperature  over  a  vertical 
height  3L  surround  the  heated  region.  This 
is  the  adjusted  temperature  at  time  t  +  ts. 

If  the  adjusted  temperature  is  rather  high,  then  the  process 

can  be  repeated  using  a  large  initial  region,  say  2L,  and 

determining  an  adjusted  temperature  at  a  later  time. 

6.6  Air-Temperature  Environment 

A  series  of  calculations  was  made  in  order  to 
evaluate  the  air- temperature  environment  and  to  determine 
the  sensitivity  of  the  environment  to  some  of  the  many 
variables  of  the  problem.  A  reference  set  of  variable 
values  was  selected  as  a  basis  for  comparison. 
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The  weapons  effect  variables  were  held  constant, 
using  a  large  weapon  detonated  at  an  optimum  height  of  burst, 
for  a  typical  overpressure  level.  The  remaining  variables 
and  their  reference  values  are  presented  in  Table  6.3. 

The  influence  of  the  heat  transfer  response  time 
is  presented  in  Figure  6.31.  It  should  be  noted  that  as 
the  response  time  decreases,  the  peak  temperature  increases 
rapidly.  This  aspect  of  the  problem  was  discussed  in  Sec¬ 
tion  6.5.3.  These  results  demonstrate  the  importance  of 
using  the  correct  value  for  this  variable.  Since  value  is 
uncertain,  it  would  appear  that  some  experimental  effort 
should  be  expended  in  order  to  reduce  the  indicated  uncer¬ 
tainty  of  the  peak  air  temperature . 

The  influence  of  a  variation  in  the  value  of  the 
transport  coefficient  is  presented  in  Figure  6.32.  Clearly 
the  reference  value  used  is  adequate  for  the  current  design 
calculations.  Figure  6.  33  illustrates  the  influence  of  the 
elevation.  The  value  of  this  design  variable  will  of  course 
be  known;  however,  it  should  be  noted  that  rather  strong 
vertical  gradients  do  exist  in  this  general  height. 

The  influence  of  the  value  of  the  land  use  factor 
is,  as  expected,  important.  This  is  demonstrated  in  Fig¬ 
ure  6.34.  Thus,  specific  site  evaluation  becomes  an  im¬ 
portant  factor  if  any  credit  for  this  effect  is  to  be 
utilized  in  the  design  of  the  dust  and  air -temperature 
conditioning  equipment. 

Another  site  variable  which  has  a  significant 
influence  on  the  air  temperature  is  the  erosion  cutoff 
depth.  Figure  6.35  demonstrates  the  variations  which  will 
occur  due  to  the  value  of  this  variable.  The  worst  case 
calculated  is  for  an  erosion  cutoff  occurring  at  a  depth  of 
0.5  inch.  This  depth  corresponds  approximately  to  the  overall 
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Table  6.3  REFERENCE  VALUES  OF  AIR-TEMPERATURE  VARIABLES 


r . . . .  . . .  1  — i 

i 

j  Variable 

i 

i 

Reference  Values 

• 

t 

'  Heat  Transfer  Response  Time  (TM) 

i 

15  ms  1 

! 

Transport  Coefficient  (K) 

0.20 

i  Erosion  Coefficient  (A) 

6  x  10° 

I  Elevation  ( Y ) 

15  ft  i 

Land  Use  Factor  (UF) 

0.5  j 

!  Erosion  Depth  Cutoff  (ZM) 

2.0  in. 

i  Bulk  Density  of  Soil  (RHOTOT) 

ioo  lb/ft J  : 

j  Density  of  Dust  Particles  (RHO) 

160  lb/ft J  | 

1  Particle  Specific  Heat  (C) 

0.20  Btu/lb  ! 

j  Soil  Conductivity  (FK) 

0.1  Btu/ft-hr-°F 

Soil  Diffusivity  (ALPHA) 

0.007  ft2/hr  ! 

Ambient  Air  Temperature  (TEMP A) 

70°F  i 

Ambient  Air  Pressure  (PRES A) 

14.7  psia  1 

— 

oos 
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Figure  6.31  Influence  of  Heat  Transfer  Response  Time 


Figure  6*32  Influence  of  Transport  Coefficient 


Figure  6.33  Influence  of  Elevation  Above  Ground 


Figure  6, 34  Influence  of  Land  Use  Factor 
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Figure  6.35  Influence  of  Erosion  Cut-off  Depth 


depth  in  which  the  thermal  energy  absorbed  by  the  soil  is 
distributed.  Additional  erosion  introduces  cold  dust  into 
the  air  and  a  significant  amount  of  cooling  occurs.  On 
the  other  hand,  if  the  cutoff  depth  is  low  (say  at  0.2  in.), 
all  of  the  thermal  energy  is  not  raised  by  the  dust  cloud 
and  a  lower  peak  temperature  thus  exists.  The  temperature 
does,  however,  remain  near  its  peak  value  for  some  time. 

A  series  of  calculations  was  made  in  which  the 
value  of  the  thermal  absorption  coefficient  was  varied. 

These  results  are  tabulated  in  Table  6.4  and  apply  to  the 
case  where  the  dust  cloud  is  (1)  very  opaque,  (2)  opaque, 
and  (3)  very  transparent.  A  value  of  0.01  lb/ft  is  used 
for  the  current  model.  For  a  very  opaque  dust  cloud,  the 
thermal  energy  incident  on  the  cloud  is  deposited  in  a  narrow 
zone  near  the  top  of  the  cloud;  hence,  a  small  temperature 
rise  occurs  early.  For  a  very  transparent  cloud,  the  thermal 
energy  is  always  deposited  near  the  ground  surface  and  the 
temperature  at  late  times  somewhat  increases.  The  intermedi¬ 
ate  case  behaves  much  like  the  very  opaque  case  where  the 
late  time  energy  absorbed  by  the  overall  cloud  is  not  felt 
at  the  lower  reaches  of  the  cloud.  In  each  case,  the 
influence  of  the  thermal  radiation  absorbed  by  the  cloud  on 
the  air  temperature  at  the  heights  of  interest  is  small. 

This  is  due  to  the  fact  that  the  rate  of  the  energy  being 
introduced  into  the  cloud  by  direct  thermal  radiation  is 
small,  compared  to  the  rate  at  which  energy  is  being  intro¬ 
duced  by  the  erosion  process.  In  the  case  illustrated, 
the  ratio  of  these  energies  is  approximately  100  to  1. 

The  influence  of  the  thermal  properties  of  the  soil 
on  the  air  temperature  was  investigated  using  the  soil 
property  values  indicated  in  Section  6.3.  These  results  are 
presented  in  Figure  6.36  and  the  curves  are  identified  by 
the  letters  a,  b,  and  c,  as  before.  The  curve  identified  by 
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the  letter  d  refers  to  the  cold  dust  case,  whereas  the  curve 
identified  by  the  letter  e  refers  to  the  no  dust  case.  The 
air  temperature  does  not  vary  greatly  for  this  normal  range 
of  thermal  properties  of  soil.  The  air  temperature  can  be 
expected  to  be  considerably  lower  for  the  case  of  very  wet 
soils.  The  influence  of  the  thermal  energy  upon  the  dust 
density  is  presented  for  these  cases  in  Figure  6.37.  The 
curve  labeled  with  the  leffer  f  corresponds  to  the  dust 
prediction  method  described  in  the  preceding  chapter.  The 
influence  of  the  heated  air  upon  the  dust  density  is  mani¬ 
fested  in  two  ways.  Fiist,  the  expansion  of  the  air  tends 
to  decrease  the  actual  dust  density  and  secondly,  the 
vertical  motion  of  the  air  lowers  the  effective  height  of 
interest,  thus  tending  to  increase  the  dust  density.  Clear¬ 
ly  the  first  process  dominates  in  the  current  situation. 

At  late  times  we  can  expect  that  the  influence  of  the  thermal 
energy  will  disappear  and  all  of  the  curves  will  approach 
the  "cold"  dust  curve. 

6.7  Multiple  Burst  Environment 

A  series  of  calculations  was  made  for  the  air- 
temperature  environment  resulting  for  specific  two-  and 

three-burst  situations.  In  each  case  the  weapons  are  of  the 
same  size  and  they  are  detonated  in  the  same  general  area. 

The  detonations  are  separated  in  time  by  approximately  one 
positive  phase  duration.  It  is  implied  that  the  respective 
fireballs  do  not  interact  with  each  other,  and  it  is  assumed 
that  the  thermal  radiation  emitted  by  the  two  or  three 
weapons,  and  incident  upon  the  site,  can  be  added  algebrai¬ 
cally.  The  yield  of  the  weapons  is  relatively  large  and 
the  site  is  located  at  an  overpressure  range  of  general 
interest . 
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The  thermal  power  curve  used  for  these  calculations 
is  presented  in  Figure  6.38.  The  dashed  line  indicates  the 
arrival  time  of  the  first  air  shock  at  the  site.  The  soil 
temperature  distribution  resulting  from  these  thermal  pulses 
is  shown  in  Figure  6.39.  The  quantity  of  energy  absorbed  by 
the  soil  decreases  fro  each  additional  weapon.  This  is  due 
to  the  successively  high  surface  temperatures  in  each  case 
and  thus,  a  greater  quantity  of  energy  is  reinitiated  back 
into  space.  The  resulting  air  temperature  curves  are  pre¬ 
sented  in  Figure  6.40  for  two  different  values  of  the  land 
use  factor.  The  height-peak  air  temperature  for  this  spe¬ 
cific  set  of  input  values  is  a  little  greater  than  900°F. 

6.8  List  of  Symbols 


erosion  coefficient, 
radiant  flux  intensity, 
specific  heat, 

slant  range  distance, 
energy , 

land  use  factor, 
buoyancy  t 

thermal  flux, 

height, 

standard  acceleration,  32.2  ft/sec^, 
heat  transfer  coefficient, 
radiation  flux, 
radiation  intensity, 

transport  coefficient, 
thermal  conductivity, 
distance, 
distance , 

number  of  particles , 
pressure , 


* 
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P ,Pe»^max 

fireball  thermal  emission, 

Qo 

thermal  flux. 

0i,0r 

thermal  power  density, 

Q»Qa»Qabs» 
Qext»^s »Qsca 

efficiencies , 

RyR^  ,Rj_ 

distance , 

Sij 

height. 

T,Te,T0 

temperature 

t  ,tg ,  t£  ,tm, 

tmax»tobs»ts 

time. 

U 

shock  velocity. 

V 

volume , 

Vi,VQ 

air  velocity 

V 

vo  lume , 

W 

weapon  yield. 

W,Wi,W0 

weight , 

w,Wij  ,wG 

weight , 

X 

distance , 

X 

distance  below  surface. 

Y 

distance , 

y 

distance , 

zc 

a 

erosion  depth, 

thermal  diffusivity,  absorption 
coefficient 

0 

coefficient  of  thermal  expansion 

c 

height , 

e 

angle,  temperature  difference. 

K 

absorption  coefficient. 

X 

wavelength. 

u 

viscosity , 

P  »P  »Pc’Psa» 

pd  density. 

a 

Stefan's  constant. 

T 

time , 

Ta 

atmospheric  transmittance . 
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XBU1 


d/d  ‘aa^oa  TBuuaqx 
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Thermal  Power  fo 


Depth,  zt  inches 

Figure  6.39  Soil  Temperature  for  Multiple  Bursts 


Temperature 


Time,  t,  seconds 


Figure  6.40  Temperature  Histories  for  Multiple  Bursts 
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Appendix  A 

COMPUTER  PROGRAM  FOR  D'JST  AND  AIR-TEMPERATURE  PREDICTION 


A.l  Discussion  of  the  Computer  Program 

The  theoretical  models  for  the  dust  and  thermal 
environment  response  to  a  simulated  nuclear  attack  condition 
have  been  coded  into  a  FORTRAN  IV  compiler  language  computer 
program  for  high-speed  numerical  prediction  calculations. 

The  program  contains  a  total  of  twelve  subprograms: 


Subprogram 
Name _ 

(Main  Program) 
WE(a-^, . . .  >a-^Q) 


_ Description 

Executive  control  of  execution,  all 
input  statements , 

Free-field  effects  for  an  air  burst, 


DUST1 

DUST2 

DUST3 

DUST4 

WEl(a^, . . .  ,ag) 
OUT  (ax) 

ATKax) 

AT2(ax) 

ST(a1,a2,a3) 


Dust  cloud  similarity  solution,  single 
burst , 

Dust  cloud  nonuniform  soil  solution, 
single  burst, 

Dust  cloud  similarity  solution,  two  or 
three  bursts, 

Late  time  dust  cloud  solution. 

Overpressure  -  Range  for  late  time 
dust  cloud, 

Output  statements  for  all  dust  cloud 
solutions , 

Air  temperature  and  dust  cloud  solution, 
uniform  ground,  single  burst, 

Air  temperature  and  dust  cloud  solution, 
biuniform  ground,  single  burst, 

Soil  temperature  solution,  and 


PE(a1,a2) 


Scaled  fireball  thermal  power  at  any  time. 


The  variables  a^,  a2,...aio  are  dummy  variables,  indicating 

i 

the  location  and  number  of  subprogram  arguments.  The  four 
distinct  dust  cloud  subroutines  and  two  distinct  air- 
temperature  subroutines  were  written  to  achieve  the  most 
economical  computational  method,  consistent  with  the  degree 
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of  information  required,  for  specific  and  general  site 
evaluation.  These  subprograms  will  be  discussed  in  detai 


A. 1.1  Executive  Main  Program 

The  main  program  is  responsible  for  the  primary 
user-oriented  flow  of  control  during  execution.  It  can  be 
regarded  as  the  executive  program  that  interfaces  a  user  s 
computational  objectives  with  the  realization  of  the  appro- 
priate  computer  results. 

All  the  input  statements  are  included  in  the  main 
program.  A  data  set  contains  values  of  all  the  input  vari¬ 
ables  necessary  for  subsequent  execution  of  one  of  the  fo  - 
lowing  basic  subroutines:  DUST1,  DUST2,  DUST3,  DUSTA,  ATI, 
or  AT2.  Any  number  of  different  data  sets  will  be  sequen¬ 
tially  executed  by  the  computer  program. 

The  first  data  card  of  any  data  set  contains  t  e 

alphanumeric  TITLE  field,  numerical  values  of  RUNT YP  and 
TOUT  and  a  logical  value  for  TAPSAV  which  are  listed  in 
Table  A.l.  The  value  of  RUNTYP  specifies  the  type  of  com- 


putation  according  to 


RUNTYP  = 


1,  DUST1  computation, 

2,  DUST2  computation, 

3,  DUST3  computation, 

4}  DUST4  computation 

5,  ATI  computation,  and 

6,  AT 2  computation. 


The  variable  IOUT  is  used  to  determine  the  level 
of  detail  of  the  computed  results  desired  as  output.  For 
example,  if  the  dust  density  is  computed  internally  at  one- 
second  intervals,  but  only  the  results  at  five-second  inter¬ 
vals  are  desired  on  the  output  listings,  then  IOUT  will  be 
five;  for  more  detailed  ouput,  the  value  of  IOUT  would  be 
suitably  decreased.  Also,  a  true  or  false  value  is  assigne 
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Table  A.l  INPUT  VARIABLES  I 


ascription  I  FORTRAN  Symbol 


Data  set  identification  (<:  72  characters)* 

TITLE 

Type  of  data  set 

RUNTYP 

Number  of  observation  time  increments 
i  between  entries  in  output  listing 

IOUT 

Binary  save  tape  option 

TAPSAV 

Number  of  observation  times 

I  MAX 

Number  of  observation  elevations 

NY 

)  Number  of  bursts* 

NB 

;  Number  of  particle  size  classes 

NCLASS 

j  Number  of  upstream  range  zones 

NZF 

i  Number  of  downstream  range  zones 

NZB 

i Time  of  second  burst  (sec)* 

TB2 

« 

;  Time  of  third  burst  (sec)* 

TB3 

Observation  time  increment  (sec) 

DIOBS 

Total  observation  time  (hr) 

TMAX 

Observation  time  increment  (hr) 

DT 

Wind  speed  (mph)* 

W  | 

I 

Wind  direction  (radians)* 

ANG 

» 

♦Included  in  output. 
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to  TAPSAV  to  indicate  whether  or  not  the  computed  results 
are  to  be  written  on  a  binary  save  tape. 

For  RUNTYP  4  (i.e.,  dust  computations  only),  the 
second  data  card  contains  values  of  other  variables  in 
Table  A.l.  The  number  of  observation  times  IMAX  must  be 
<;  200.  The  number  of  elevations  NY  to  be  entered  as  input 
data  must  be  £  3.  The  dust  density  history  is  computed  at 
each  of  the  specified  elevations. 

For  RUNTYP  «  1,  2,  or  4  the  value  of  NB  must  be 
1,  i.e.,  single-burst  dust  computations.  In  these  cases, 
the  time-of-burst  variables  TB2  and  TB3  are  superfluous; 
but,  for  the  purpose  of  standardized  procedure,  any  such 
variables  that  may  have  assumed  arbitrary  values  will  be 
set  equal  to  zero.  For  NTYP  =  3,  NB  must  be  2  or  3,  giving 
the  number  of  bursts  for  the  multiple -burst  dust  computation; 
if  NB  =  2,  then  TB2  must  be  specified  and  TB3  =  0;  and  if 
NB  =  3,  then  both  TB2  and  TB3  must  be  specified. 

The  number  of  particle  classes  NCIASS  must  be  <:  8. 
For  NTYP  =  1,  3,  or  4  the  range  zone  numbers  NZF  and  NZB 
must  be  zero.  For  NTYP  =  2,  the  values  of  NZF  and  NZB  must 
be  specified  and  less  than  21. 

The  value  of  the  observation  time  increment  DTOBS 
is  arbitrary;  the  product  of  IMAX  and  DTOBS  gives  the  time 
corresponding  to  the  end  of  computations. 

For  RUNTYP  =  4,  the  remaining  variables  in  Table  A.l 
are  provided  on  the  next  data  card.  These  are  variables  that 
relate  only  to  the  late  time  dust  computations  and  include 
TMAX,  the  time  corresponding  to  the  end  of  computations,  DT, 
the  time  step,  W,  the  wind  speed,  and  ANG,  the  wind  direction. 
Note  that  for  RUNTYP  =  4  the  variables  IMAX  and  DTOBS  are 
superfluous. 


For  RUNTYP  =  5  or  6,  i.e.,  air  temperature  compu¬ 
tations,  the  second  data  card  contains  values  of  the  vari¬ 
ables  given  in  Table  A. 2.  The  variables  UF  and  ZM  pertain 
only  to  the  uniform  ground  solution  and,  therefore,  have 
values  only  if  RUNTYP  =  5.  Similarly,  the  variables  PI, 

P2,  and  RC  pertain  to  the  biuniform  ground  solution  and 
have  values  only  if  RUNTYP  =6. 

The  next  data  card  of  each  set  contains  values 
for  the  soil  particle  density,  bulk  soil  density,  ambient 
air  temperature  and  pressure,  and  the  transport  and  erosion 
constants  which  are  given  in  Table  A. 3.  The  other  variables 
in  this  table  are  discussed  below. 

Five  variables  establish  the  nuclear  attack  con¬ 
dition;  the  weapon  yield,  the  overpressure  at  the  plant  site, 
and  three  height -of -burst  variables  that  select  among  optional 
choices  for  the  height  of  burst.  For  RUNTYP  f  3,  the  only 
burst  variables  entered  as  input  data  are  YIELD,  PRESO,  HOB, 
HOBS,  and  HOBO.  For  RUNTYP  «  5  or  6,  this  will  be  the  last 
data  card  of  the  data  set.  For  RUNTYP  -  3,  the  burst  vari¬ 
ables  entered  as  input  data  are  YIELD 1,  PRES01,  H0B1,  H0BS1, 
and  H0B01  —  and  YIELD2 ,  PRES02 ,  H0B2,  H0BS2 ,  and  H0B02  for 
NB  =  2  or  3,  and  in  addition,  YIELD3 ,  PRES03 ,  H0B3,  HOBS 3, 
and  H0B03  for  NB  =  3. 

The  data  card  following  the  burst  variables  card(s) 
contains  the  one-dimensional  array  for  the  particle  diameters 
D(N) ,  N  *  1,  NCLA.SS ,  for  RUNTYP  <;  4. 

The  variable  NZR1,  which  equals  the  number  of  up¬ 
stream  zones  plus  the  number  of  downstream  zones  plus  one 
(the  zone  containing  the  site),  is  then  computed.  The  next 
NZR1  data  cards  contain  the  two-dimensional  array  for  the 
adjusted  weight  fractions,  P(N,NN),  N  =  1,NCLASS,  NN  *  1,NZR1, 
where  each  card  contains  NCIASS  values  of  P  for  the  respective 
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Table  A.’  INPUT  VARIABLES  II 


Description 


FORTRAN  Symbol 


Soil  conductivity  (Btu/ft-hr-°F)* 

Soil  diffusivity  (ft^/hr)* 

Observation  time  factor 

Soil  particle  specific  heat  (Btu/lb)* 

Heat  transfer  response  time  (msec)* 
Observation  elevation  (ft)* 

Uniform  ground  land  use  factor* 

Uniform  ground  erosion  depth  cutoff  (in.)* 
Biuniform  ground  near  land  use  factor* 

•  Biuniform  ground  far  land  use  factor* 
Biuniform  ground  range  boundary  (ft)* 


ALPHA 

BETA 


*Included  in  output. 


Table  A. 3  INPUT  VARIABLES  III 


■  .  : — -  —  _  .  —  ■  ■  ■■  i 

Description 

FORTRAN  Symbol 

3  * 

Density  cf  soil  particles  (lb/ft  ) 

RHO 

Bulk  density  of  soil  (lb/ft  ; 

RHOTOT 

Ambient  air  temperature  (°F) 

TEMPA 

Ambient  air  pressure  (psia) 

PRESA 

Transport  constant 

KH 

Erosion  constant 

A 

■ 

Single  Burst 

Weapon  yield  (MT) 

YIELD 

Overpressure  (psig) 

PRESO 

*  Height  of  burst  (ft) 

HOB 

Scaled  heif\.  of  burst  (ft /Mr'  ) 

HOBS 

Optimum  height  of  burst 

HOBO  j 

l 

Multiple  Bursts 

i 

Weapon  yield  (MT) 

YIELDl,  YIELD 2,  YIELD3 

Overpressures  (psi) 

PRES01,  PRESO 2 ,  PRES03  ! 

Heights  of  burst  (ft) 

H0B1,  H0B2,  HOB 3 

!  Scaled  heights  of  burst  (ft /MI1'; 

H0BS1,  HOBS 2,  HOBS 3 

j  Optimum  heights  of  burst 

H0B01,  H0B02,  H0B03 

| 

Particle  diameters  (mm) 

D(N) ,  N=1,NCLASS 

Adjusted  weight  fractions 

P(N,NN),  N=1,NCLASS, 

NN=1,NZR1 

k 

Observation  elevations  (ft) 

Y(N),  N=1,NY 

Forward  zone  boundaries  (ft) 

R.ZF  (N) ,  N=1,NZF 

Backward  zone  boundaries  (ft) 

_ 

RZB(N),  N=1,NZB 

i - - - — —  : ~ — —  - .  i 

“fc 

Included  in  output. 

> 
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zone  NN.  for  RUNTYF  =  1,  3  or  4,  a  uniform  ground  condition 
is  assumed,  and  therefore  NZR1  =  1;  for  RUNTYP  =  2,  NZR1  <  41. 

For  RUNTYP  <  4,  the  specified  observation  elevations 
Y(N),  N  =  1,NY  are  entered  from  the  next  data  card  that  can, 
at  most,  specify  three  elevations.  If  RUNTYP  =  1,  3,  or  4, 
this  is  the  last  data  card  of  the  data  set.  If,  however, 
RUNTYP  =  2,  then  the  forward  and  backward  zone  ooundaries 
RZF(N),  N  =  1,NZF  and  RZB(N),  N  =  1,1  .B  are  entered  from  the 
next  two  and  final  data  cards,  respectively.  The  weight 
fractions  and  zone  boundaries  are  introduced  in  the  order  of 
increasing  distance  in  front  of,  and  then  in  the  order  of 
increasing  distance  behind,  the  plant  location.  The  plant 
location  is  not  considered  a  zone  boundary. 

Sample  problem  data  sets  for  a  sample  site  are 
presented  in  Section  A. 2. 3  They  provide  examples  of  the 
data  card  ordering  proced*  discussed  above. 

The  main  progi  then  directs  control  to  the  compu¬ 
tational  subroutines  depending  on  the  value  of  RUNTYP  accord¬ 
ing  to 


RUNTYP  = 

i 


1,  CALL  DUST1, 

2,  CALL  DUST2, 

3,  CALL  DU ST 3 , 

4,  CALL  DU ST 4 , 

5,  CALL  ATI,  and 


6 ,  CALL  AT 2 

For  RUNTYP  =  1,  2,  or  3,  upon  return  of  control 
from  the  appropriate  subroutine,  the  output  subroutine  OUT 
is  called;  this  provides  the  off-line  listing  of  computa¬ 
tional  dust  cloud  results.  The  TAFSA.V  variable  is  then 
tested  and,  if  it  is  true,  a  binary  save  tape  of  these 
results  is  generated.  For  RUNTYP  =4,  5,  or  6 ,  the  off¬ 
line  listing  is  generated  in  either  the  DUST4,  ATI  or 
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AT2  subprogram,  respectively;  the  TAPSAV  option  is  similarly 
tested  for  binary  storage  of  the  temperature  computations. 

A. 1.2  Weapons  Effects  Subprogram 

Th^  dust  and  thermal  environment  prediction  program 
requires  certain  weapons  effects  variable  values  as  data. 

These  values,  in  turn,  must  be  determined  from  a  given  nuclear 
attack  condition.  A  weapon  effects  subroutine  was  written  to 
provide  these  variable  values  over  the  range  of  variables  of 
interest . 

The  many  weapon  effects  variables  are  interrelated 
such  that  the  attack  situation  can  be  uniquely  defined  by 
specifying  a  limited  number  of  variable  values.  To  limit  the 
complexity  of  the  weapons  effects  subroutine,  the  defining 
weapon  effects  variables  were  restricted  to  weapon  yield, 
YIELD(MT),  overpressure,  PRESO  (psi),  and  height  of  burst  or 
equivalent. 

The  height  of  burst  option  permits  us  to  choose 
(1)  the  actual  height  of  burst  HOB  (ft),  (2)  the  scaled  height 
of  burst  at  standard  air  conditions  HOBS  (ft/MT1^),  (3)  the 
optimum  height  of  burst  HOBO,  or  (4)  an  overhead  burst.  The 
subroutine  sequentially  examines  the  value  of  HOB,  HOBS,  and 
HOBO  until  an  acceptable  value  is  indicated.  A  negative  HOB 
value  is  not  acceptable;  hence  an  input  value  such  as  -1  would 
cause  the  subroutine  to  bypass  this  option.  Any  positive  value, 
including  zero,  is  an  acceptable  value.  A  negative  or  zero 
value  for  HOBS  is  not  acceptable,  and  a  value  less  than  0.1 
is  not  acceptable  for  HOBO.  If  these  three  input  values  are 
not  acceptable,  then  option  (4)  is  automatically  selected. 

The  final  input  data  needed  to  define  the  attack  situation 
are  the  ambient  pressure  PRESA  and  the  temperature  TEMPA. 

This  subroutine  is  entered  through  a  CALL  WE(a^,...,a^Q) 
statement  in  the  main  program  where  a^  represents  the  arguments 
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of  the  subroutine.  Execution  then  proceeds  through  four  parts: 
(1)  an  initial  computation  in  which  a  variety  of  scale  factors 
and  other  numbers  are  generated  for  subsequent  use,  (2)  the 
calculation  of  those  weapon  effects  variables  which  are  not  a 
function  of  the  height  of  burst,  (3)  the  calculation  of  those 
variables  which  are  a  function  of  the  height  of  burst ,  and 

(4)  the  calculation  of  some  thermal  data.  The  first  part 
consists  of  six  steps: 

1.  Compute  absolute  ambient  temperature,  TEMP  AS. 

2.  Compute  scale  factors  for  yield,  CRY,  ambient  pressure, 

PSF,  ambient  temperature ,  TSF,  and  overpressure,  SOP. 

3.  Compute  ambient  sound  velocity,  CVELA. 

4.  Compute  shock  strength,  XI. 

5.  Compute  two  height -of -burst  scale  factors,  SRB  and  SRC. 

6.  Compute  square  root  of  yield,  SQRTW. 

The  second  part  of  the  subroutine  consists  of  seven  steps: 

1.  Compute  time  of  maximum  thermal  emission,  TMAX. 

2.  Compute  exponent  for  velocity  pulse,  KV.  If  KV  is  less 
than  1,  set  KV  equal  to  1. 

3.  Compute  positive  phase  duration  of  velocity  pulse,  TO, 

•  using  one  of  two  formulas  depending  upon  the  value  of  SOP. 

4.  Compute' the  air  temperature  at  the  shock  front,  TES. 

5.  Compute  the  air  density  at  the  shock  front,  RHOS,  and  an 
average  value  within  the  blast  wave,  RHOF. 

6.  Compute  the  positive  phase  duration  of  the  overpressure, 
TOP. 

7.  Compute  three  coefficients  for  the  overpressure  waveform, 
CPI,  CP2,  and  CP3 . 

The  third  part  consists  of  four  main  steps  which  are  subdivided 
as  follows: 

1.  This  step  examines  the  height -of -burst  input  information 
and  determines  whether  the  height  of  burst  is  such  that 
the  following  shock  regimes  exist: 

a.  Mach  reflection, 

b.  regular  reflection, 

c.  overhead  burst. 
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It  also  computes  the  actual  height  of  burst,  HOB.  or  the 
scaled  height  of  bu:*:t _  SHOB,  whichever  is  not  given  as 
input . 

2.  This  step  computes  the  weapons  effects  variables  if  the  Mach 
reflection  regime  applies.  The  variables  which  are  evalu¬ 
ated  are: 

a.  shock  velocity,  U, 

b.  peak  air  velocity,  VO, 

c.  time  of  arrival  of  shock  at  site,  TASITE, 

d.  scaled  range,  SRA,  and 

e.  angle  of  incidence,  THETA  (if  theta  is  very  small, 
set  theta  equal  to  rr/2. 

3.  This  step  computes  the  weapons  effects  variables  if  the 
regular  reflection  regime  applies.  The  following  calcu¬ 
lation  or  tests  are  made: 

a.  compute  maximum  permissible  height  of  burst,  SRD, 

b.  compare  SHOB  with  SRD;  if  SHOB  is  greater  than  or 
equal  to  SRD,  transfer  to  overhead  burst  calcula¬ 
tions  (this  automatically  sets  SHOB  to  a  permis¬ 
sible  value) , 

c.  compute  factor,  CF,  for  velocity  calculations, 

d.  compute  angle  of  incidence,  THETA, 

e.  compute  shock  velocity,  U, 

f.  compute  peak  air  velocity,  VO, 

g.  compute  scaled  height  of  burst,  SRE, 

h.  compare  SHOB  with  SRE;  use  one  of  two  equations  to 
compute  time  of  arrival,  TASITE. 

4.  This  step  computes  the  weapons  effects  variables  for  an 
overhead  burst:  shock  velocity,  U,  peak  air  velocity, 

VO,  angle  of  incidence,  THETA,  height  of  burst,  HOB,  and 
time  of  arrival,  TASITE. 

5.  The  last  step: 

a.  computes  cosine  of  theta,  COSTH, 

b.  computes  slant  range,  SRANGE, 

c.  examines  magnitude  of  HOB  (if  HOB  is  small,  computes 
slant  range  from  a  second  equation), 

d.  computes  an  exponential  factor,  OPT,  from  tabular 
data  after  first  determining  the  height-of-burst 
interval  applying  to  the  current  case, 

e.  computes  thermal  attenuation  factor,  TAU,  and 

f.  computes  peak  radiant  power  incident  at  site,  CONSV>T1.# 

Control  is  then  returned  to  the  main  program  for  subsequent 
execution. 

A. 1.3  Uniform  Ground  Dust  Subprogram 

This  single-burst  uniform-ground-condition  subroutine 
makes  use  of  the  similarity  solution  as  well  as  the  analytical 
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expression  tor  cloud  height.  Thus,  the  ceding  reduces  to  a 
simple  evaluation  of  these  equations  for  a  set  of  transport 
times.  The  dust  densities  for  various  particle  size  classes 
and  for  various  elevations  are  computed  at  uniform  time  inter¬ 
vals  measured  at  a  fixed  range  (the  plant  location)  to  which 
the  corresponding  transport  times  are  related. 

This  subroutine  is  entered  through  the  CALL  DUST1 
statement  in  the  main  program,  the  calling  program.  Execution 
then  proceeds  through  two  main  computational  parts  of  this 
subroutine:  (1)  an  initial  computation  phase  in  which  a 

variety  of  numbers  is  generated  for  subsequent  use,  and  (2)  an 
evaluation  of  the  cloud  height  and  the  computation  of  the  dust 
density  from  the  similarity  solution.  The  first  part  consists 
of  four  steps: 

1.  Computation  of  the  maximum  air  displacement  RO  and  the 
related  time  factors,  transit  time  TAO,  modified  positive- 
phase  duration  TOM,  and  the  exponential  factor  TAU. 

2.  Computation  of  terminal  velocities  for  each  particle  size 
class  VT(N) . 

3.  If  VO  is  less  than  10  fps,  CUTOFF  is  set  equal  to  false 
and  control  returns  to  the  main  program,  thereby  elimi¬ 
nating  any  further  computation  of  the  specific  problem. 

4.  Evaluation  of  a  number  of  constants. 

The  second  part  consists  of  eight  steps: 

1.  Initialization  of  all  dust  density  storage  locations  to 
zero. 

2.  Computation  of  observation  times  TOBS(I). 

3.  Computation  by  iteration,  of  the  air  displacement  R 
corresponding  to  TOBS(I).  Evaluate  the  corresponding 
transit  time  T. 

4.  Computation  of  the  height  of  the  cloud  HST  in  the  absence 
of  gravity. 

5.  Computation  of  the  height:  of  the  cloud  H  for  each  particle 
class  (and  for  each  time). 

6.  Computation  of  the  value  of  the  similarity  variable  X  for 
each  elevation  Y(L)  and  for  each  particle  size  class  and 
for  each  time. 
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/.  Evaluation  of  the  dust  density  contribution  DYC(I,N,L) 
from  the  similarity  solution. 

8.  Summation  over  all  particle  size  classes  to  obtain  the 
total  dust  density  D(I,L)  at  each  elevation  for  each 
TOBS(I). 

Control  is  then  returned  to  the  main  program  for 
subsequent  output  of  the  computed  variables  and  execution 
of  the  remaining  problem  data  sets. 

A. 1.4  Variable  Ground-Dust  Subprogram 

The  single-burst  variable -ground  codition  sub¬ 
routine  computes  the  dust  density  variation  at  various 
height  and  times  in  accordance  with  the  basic  IITRI  dust 
model.  Equation  5.13  is  used  to  integrate  all  of  the 
elemental  dust  contributions  picked  up  along  that  particular 
air  column  which  arrives  at  the  site  location  at  the  time 
of  interest.  The  total  dust  density  DY  and  the  dust  density 
for  each  particle  size  class  DYC  are  computed  and  printed 
out  for  each  elevation  Y. 

This  subroutine  is  entered  through  the  CALL  DUST2 
statement  in  the  main  program.  Execution  then  proceeds 
through  two  main  computational  parts:  (1)  an  initial  compu¬ 
tation  in  which  a  variety  of  numbers  are  generated  for  sub¬ 
sequent  use,  (2)  the  integration  or  summation  in  which  the 
elementary  dust  density  contributions  are  accumulated  for 
each  increment  of  time  and  along  each  integration  path, 
and  for  each  particle  size  class.  The  initial  part  consists 
of  12  steps: 

1.  Computation  of  the  maximum  air  displacement  RO  and  the 
related  time  factors,  transit  time  TAO,  the  modified 
positive-phase  duration  TOM,  and  the  exponential  factor 
TAU . 

2.  Computation  of  terminal  velocities  for  each  particle 
size  class  VT(N). 

3.  If  VO  is  less  than  10  fps,  CUTOFF  is  set  equal  to  false 
and  control  returns  to  the  main  program,  thereby  elimi¬ 
nating  any  further  computation  of  the  specific  problem. 
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4.  Computation  of  the  time  increment  DT  for  integration 
along  each  particle  path. 

5.  Evalutation  of  constants. 

6.  Conversion  of  range  boundaries  RZF  and  RZB  into  an 
ordered  set  of  range  constants  R2(I)  measured  from 
R  =  0  (the  site  is  located  at  R  =  RO) . 

7.  Computation  of  the  observation  time  TOBS(I)  and  a 
corresponding  constant  CONS. 

b.  Computation,  by  iteration,  of  the  range  R1(I)  from  which 
the  air  column  started  and  just  reached  the  site  at  the 
corresponding  observation  time. 

9.  Computation  of  the  shock  arrival  time  at  R1(I),  the  total 
transit  time  DT1(I),  and  an  exponential  factor  TAU1. 

10.  Computation  from  one  of  two  equations,  of  the  height  of 
the  dust  cloud  H1(I)  in  the  absence  of  gravity  at  the 
corresponding  observation  time. 

11.  Computation,  for  each  of  500  integration  intervals,  of 
the  transit  time  TI(J);  evaluate  two  corresponding 
constants;  and  then  computation  of  the  velocity  VEL(J) 
and  the  air  displacement  DIS(J). 

12.  Computation,  from  one  of  two  equations,  of  the  height 
of  the  cloud  HSTAR(J)  for  each  integration  element. 

The  second  part  of  the  program  consists  of  four  steps: 

1.  Initialization  of  all  dust  density  storage  locations 
to  zero. 

2.  For  each  integration  path  (i=l,IMAX),  for  each  particle 
size  class  (M=I,NCLASS)  and  for  each  element  (J=1,JMAX) 
along  the  Ith  integration  path,  computation  of 

a.  the  erosion  factor  FACT:  compare  with  zero,  and 
set  equal  to  zero  or  unity; 

b.  the  range  R  of  Jth  element;  compare  with  range 
constants  R2(K),  determine  ordered  zone  number  K 
and  convert  to  zone  number  L,  which  corresponds 
to  order  of  input  data  P(M,L); 

c.  the  increment  of  mass  of  the  dust  packet  DM; 

d.  the  time  increment  after  dust  pickup  DELT ;  and 

e.  the  height  of  the  dust  cloud  H,  corresponding  to 
the  given  packet  of  dust. 

3.  For  each  elevation  Y(L),  computation  of  the  incremental 
dust  contributions  and  integration  over  all  J  elements 
to  obtain  DYC(I,M,L). 

4.  Summation  of  the  dust  contribution  from  each  particle 
size  class  to  obtain  the  total  dust  density  DY(I,L)  at 
each  elevation  Y(L)  and  for  each  observation  time 
TOBS(I)  . 
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The  range  boundaries  are  then  recomputed  and  the 
adjusted  weight  fractions  reindexed  for  output  format  pur¬ 
poses.  Control  is  then  returned  to  the  main  program. 

A. 1.5  DUST3  Multiple  Burst-Dust  Subprogram 

The  multiburst  subroutine  was  written  for  a  uniform 
ground  condition  and  represents  a  first  estimate  of  the 
influence  that  a  multiweapon  attack  has  upon  the  dust  envi¬ 
ronment,  which  is  created  by  the  resulting  airblast  winds. 

The  multiweapon  attack  problem  introduces  a  large  number  of 
additional  weapon,  timing,  and  geometrical  variables  into  an 
already  complex  problem.  To  make  the  problem  manageable  at 
this  stage  of  development  of  the  dust-environn ^nt  predictions, 
a  number  of  simplifying  assumptions  were  made.  The  assumption 
or  use  of  a  uniform  ground  condition  represents  one  such 
simplification.  This  assumption  enables  us  to  disregard  the 
settling  or  fallout  of  the  dust  from  the  first  burst  and  its 
special  redistribution  over  the  general  area  of  the  site.  It 
is  also  consistent  with  the  basic  assumption  that  there  is  no 
interaction  between  particles  of  different  size  classes  and 
their  respective  availability.  Thus,  the  layering  of  soil 
of  one  distribution  on  the  soil  of  another  distribution  need 
not  be  considered. 

The  interaction  of  one  air  blast  field  with  another 
represents  a  very  complex  nonlinear  interaction  problem,  ex¬ 
cept  wher,  the  intensity  of  one  of  the  fields  (the  first  one) 
lias  decayed  somewhat.  Thus,  it  is  essential  to  require  that 
the  time  interval  at  the  site  location  between  the  arrival  of 
the  respective  shock  parts  be  long  enough  to  allow  the  air  to 
approach  a  near-quiescent  state.  Such  a  state  first  occurs 
at  the  time  of  the  positive-phase  duration  of  the  overpressure 
and  exists  thereafter  for  that  particular  burst. 
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Under  this  condition,  the  velocity  fields  generated 
by  the  various  bursts  can  be  added  vectorially.  To  simplify 
the  computational  procedure  further,  it  was  assumed  that  the 
various  bursts  occur  in  the  same  general  direction  from  the 
site.  Hence,  the  velocity  fields  can  be  now  treated  as  scalar 
quantities  and  added  algebraically.  It  is  obvious  that  the 
addition  of  the  positive  phase  of  the  second  burst  to  the 
negative  phase  of  the  first  burst  will  reduce  the  air  velocity 
and,  hence,  reduce  the  dust  density.  This  degradation  in  air 
velocity  should  be  small  relative  to  the  primary  effect  of 
the  second  burst.  If  the  two  weapons  are  detonated  in  oppo¬ 
site  directions  relative  to  the  site,  the  velocity  would  be 
somewhat  enhanced;  however,  such  an  attack  environment  is 
not  considered  very  probable. 

The  multiburst  subroutine  for  uniform  soil  con¬ 
ditions  is  very  similar  to  th*  corresponding  single-burst 
subroutine  in  that  use  is  made  of  the  similarity  solution. 
However,  the  cloud  height  in  the  absence  of  gravity  effect 
is  computed  by  direct  numerical  integration  of  the  control¬ 
ling  differential  equation,  the  first  term  of  Equation  5.4. 

A  classical  three-point  integration  formula  was  used  within 
a  given  time  interval.  The  multiburst  subroutine  treats 
either  two  or  three  bursts.  The  shock  velocities  U  for  the 
various  bursts  are  averaged  for  computational  purposes  as 
are  the  air  densities  RHOF.  This  subroutine  is  entered 
through  the  CALL  DFST3  statement  in  the  main  program. 

Execution  cnen  proceeds  through  two  computational  parts, 
the  first  dealing  with  the  evaluation  of  a  variety  ol 
numbers,  and  the  second  dealing  with  the  density  computation 
from  the  similarity  solution.  A  series  of  logical  constants 
is  used  to  determine  which  of  three  possible  equations  will 
be  used  to  compute  the  cloud  height.  This  is  necessary  to 
avoid  the  effect  of  negative  cloud  heights.  The  first  part 
consists  of  13  steps: 
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1.  Definition  of  arithmetic  statement  functions  for  the 
velocity  and  displacement  contribution  from  the  varic is 
bursts . 

2.  Computation  of  average  values  of  U  and  RHOF. 

3.  Computation  of  terminal  velocity  for  each  particle  size 
class  VT (N) . 

4.  Computation  of  the  maximum  air  displacements  R01  and 
R02  for  the  first  two  bursts . 

5.  Computation  of  the  corresponding  time  factors  T0M1, 

T0M2,  and  a  number  of  constants  FK11  to  FK42. 

6.  Computation  of  constant  C0NS1,  accuracy  criteria  CRIT1, 
and  maximum  observation  tine  TMAX. 

7.  Computation  of  the  observation  times  TOBS(I). 

8.  Computation,  by  iteration,  of  the  transit  time  TCI 
between  the  start  of  the  flow  and  the  time  the  second 
shock  point  reaches  the  air  column. 

9.  If  TCI  >  TMAX,  CUTOFF  is  set  equal  to  false,  and  control 
returns  to  the  main  program,  thereby  terminating  this 
problem. 

10.  Evaluation  of  critical  values  of  I,  IB2,  IB21,  and  IB22, 
separating  the  first  and  second  bursts. 

11.  When  three  bursts  are  used,  computation  of 

a.  R03 ,  T0M3 ,  FK13,  FK23,  and  FK43, 

b.  the  time  of  maximum  displacement  for  the  second 
burst  by  iteration, 

c.  the  transit  time  TC2  by  iteration,  and 

d.  if  TC2  >  TMAX,  CUTOFF  is  set  equal  to  false  and 
control  return -  to  the  main  program,  and 

e.  critical  values  of  I(IB3,  IB31,  and  IB32). 

12.  Computation  of  transit  times  TA(I)  from  one  of  three 
equations,  depending  upon  the  burst  interval.  When 
only  two  bursts  occur,  IB3  is  set  equal  to  IM4X. 

13.  Computation  of  the  height  of  the  cloud  HSTAR(I)  in  the 
absence  of  gravity  for  each  transit  time  TA(I).  This 
is  done  over  six  regions  for  1=1,  1-2,  to  IB2,  I=IB2, 
I=IB21  to  IB3 ,  I=IB31,  and  I=IB32  to  IMAX.  Critical 
values  of  H,  (HC1  and  HC2)  are  also  computed,  corre¬ 
sponding  to  TCI  and  TC2 . 

The  second  portion  of  the  program  consists  of  four  steps: 

1.  Initialization  of  all  dust-density  storage  locations  to 
zero. 

2.  Initialization  of  logical  constants  Ql,  Q2 ,  and  Q3  to 
TRUE . 
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For  each  particle  size  class, 

a.  compute  the  cloud  height  H  for  each  time  TOBS(I). 

b.  compare  H  with  zero. 

c.  if  H  is  greater  than  zero,  compute  value  of 
similarity  variable  and  then  compute  density 
contribution  DCY(I,N,L)  for  each  elevation. 

d.  if  H  ^  0,  set  Q1  *  FALSE,  set  I  to  IB21. 

e.  compute  the  cloud  height  H  for  each  time  TOBS(I), 
I=IB21  to  IB3  from  one  of  two  equations,  depend¬ 
ing  upon  logical  constant  Ql. 

f.  compare  H  with  zero. 

g.  if  H  is  greater  than  zero,  compute -value  of 
similarity  variable  and  then  compute  the  density 
contribution  DCY(I,N,L)  for  each  elevation. 

h.  if  H  <  o,  set  Q2  =  FALSE,  set  I  to  IB31  (or  stops 
for  two  bursts). 

i.  compute  the  cloud  height  H  for  each  time  TOBS(I), 
I=IB31  to  IM\X  from  one  of  three  equations, 
depending  upon  the  logical  constants  Ql  and  Q2. 

j .  compare  H  with  zero . 

k.  if  H  is  greater  than  zero,  compute  value  of 
similarity  variable,  then  compute  the  density 
contribution  DCY(I,N,L)  for  each  elevation. 

l.  if  H  <  0,  set  Q3  =  FALSE.  This  terminates  the 
calculations  for  the  Nth  particle  size  class. 

(\.  Summation  of  the  dust  contribution  for  each  particle  size 
class  to  obtain  the  total  density  DY(J,L)  at  each  eleva¬ 
tion  Y(L)  and  for  each  observation  time  TOBS(I). 

Control  is  then  returned  to  the  main  program. 

A. 1.6  Late  Time  Dust  Subprogram 

This  single-burst  uniform-ground-condition  sub¬ 
routine  makes  use  of  two  modified  forms  of  the  simulation 
solution  to  compute  the  dust  density  at  late  times.  One 
solution  relates  to  the  growth  phase  of  cloud  development 
whereas  the  other  solution  relates  to  the  settling  phase. 

The  coding  consists  of  determining  the  apparent  location 
of  the  site  within  the  overall  dust  cloud  formation  at 
uniform  time  intervals  and  then  evaluating  the  corresponding 
overpressure  levels,  the  blast  wave  variables,  and  finally, 
the  maximum  cloud  height  at  these  locations.  The  dust 
densities  for  various  particle  size  classes  are  then 


A-21 


evaluated  by  determining  which  cloud  development  phase 
cpplies  and  using  the  appropriate  modified  simularity 
solution.  Due  to  limitations  of  the  weapons  effects  equa¬ 
tions  the  dust  densities  at  locations  for  which  the  over¬ 
pressure  exceeds  60  psig  are  set  equal  to  a  very  large 
number  and  therefore  printed  out  as  "******". 

The  subroutine  is  entered  through  a  CALL  DUST4 
statement  in  the  main  program.  A  series  of  four  function 
statements  is  defined  in  this  subroutine  which  relate  to 
the  coordinates  of  constant  overpressure  lines  in  the 
height  of  burst-range  space.  Execution  then  proceeds 
through  four  main  computational  parts:  (1)  a  preliminary 
calculation  part  in  which  certain  numbers  are  generated 
for  subsequent  use,  (2)  an  initialization  part  in  which 
the  individual  dust  density  contributions  are  set  equal 
to  zero,  (3)  an  iteration  procedure  in  which  a  critical 
range  is  evaluated,  and  (4)  the  basic  dust  density  evalu¬ 
ations  which  are  made  for  each  time  step.  The  initial 
part  consists  of  the  following  steps: 

1.  Set  the  value  of  IC  and  Q,  both  of  which  are  used  to 
control  the  computational  flow  when  the  excessive 
overpressure  override  is  needed. 

2.  Compute  the  number  of  computational  or  time  steps  IM. 
Compare  with  the  allowable  dimension  size  and  limit 
to  201  if  exceeded. 

3.  Compute  the  fluid  density  RHOF  and  the  terminal 
velocity  VT(N)  for  each  of  the  particle  size  classes. 

4.  Evaluate  a  series  of  constants,  CRY,  PSF,  and  TSF. 
Compute  the  scaled  overpressure  S. 

5.  Compare  HOB  or  other  height  of  burst  input  data  with 
appropriate  criteria  to  determine  which  reflection 
region  is  applicable  for  the  starting  calculation. 

Then  compute  the  range  of  the  site,  RS. 

6.  Compute  a  series  of  constants,  C0NS1  and  RR.  Set  the 
value  of  the  overpressure  PR  equal  to  PRESO. 

The  second  part  consists  of  equating  the  total  dust  densi¬ 
ties  DY(I,L)  and  the  particle  class  contribution  DYC(I,N,6) 
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to  zero  for  each  height,  L,  of  interest  and  each  time,  I. 

The  third  part  consists  of  the  following  steps: 

1.  Set  the  scaled  overpressure  SP  equal  to  15. 

2.  Using  a  Newton-Ralphson  iteration  procedure  solve  for 
which  the  actual  height  of  burst  is  the  optimum  height 
of  burst.  An  accuracy  criteria  of  0.05  psi  is  used  on 
the  convergence  interval. 

3.  Compare  SP  with  0.1  and  if  larger  use  function  RCF(SP) 
to  evaluate  critical  range  RCRIT.  Otherwise  set 
RCRIT  =  1000000. 

The  final  part  consists  of  the  following  steps: 

1.  Compute  the  observation  time  TOBS(I),  the  travel  dis¬ 
tance  RR  and  the  range  R  from  the  dust  cloud  origin. 
Compute  the  initial  value  of  the  scaled  overpressure 
S  based  upon  the  last  overpressure  evaluation. 

2.  Test  IC  and  Q  to  determine  if  overpressure  override 
applies . 

3.  Compare  R  with  RCRIT  and  select  one  of  two  overpressure 
functions  with  wl  ich  to  use  with  a  Newton-Ralphson 
iteration  procedure.  An  accuracy  criterion  of  0.1  psi 
is  used  on  the  convergence  interval.  Evaluate  the 
overpressure  PR. 

4.  Compare  PR  with  60  psi  to  determine  if  the  overpressure 
override  should  be  activated.  If  so  go  to  statement  121. 

5.  CALL  weapons  effects  subprogram  WEI  (PR)  to  evaluate  maxi¬ 
mum  cloud  height  HST. 

6.  Compute  change  in  cloud  height  DH. 

7 .  Compare  DH  with  0  to  determine  which  of  two  modified 
simularity  solutions  are  to  be  used  to  evaluate  the  dust 
density  contributions.  Compare  height  of  cloud  with 
elevation  Y(L).  If  elevation  is  larger  than  bypass 
calculation,  otherwise  evaluate  similarity  variable  X. 
Compute  dust  density  contribution  D¥C(I,N,L). 

8.  Statement  121  activates  the  overpressure  override  by 
setting  Q  =  FALSE  and  evaluating  the  critical  index  IC 
where  the  override  is  to  be  removed. 

9.  This  group  of  statements  are  used  to  evaluate  the  dust 
density  contribution  D5fC(I,N,L)  where  the  overpressure 
override  condition  is  active. 

10.  The  total  dust  densities  DY(I,L)  are  obtained  by  summing 
over  the  particle  size  classes  for  each  elevation  Y(L) 
and  observation  time  TOBS(I). 

Control  is  then  returned  to  the  main  program. 
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A. 1.7 


WEI  (Late  Time  Dast  Weapons  Effects)  Subprogram 

The  late  time  dust  subprogram  require  that  certain 
blast  wave  variable  be  evaluated  so  that  the  maximum  cloud 
height  can  be  established  for  each  time  step.  The  variables 
are  all  functions  of  the  overpressure  level  which  is  evalu¬ 
ated  in  the  subprogram.  This  subroutine  is,  therefore, 
entered  through  a  CALL  WEI  (a^,...a9)  statement  in  the  Late 
Time  Dust  Subprogram.  It  consists  of  the  following  steps: 

1.  Compute  the  scaled  overpressure  SO  and  shock  strength  XI. 

2.  Compute  the  exponential  coefficient  KV.  Compare  with 
1  and  if  less  set  KV  =  1. 

3.  Compute  positive  phase  duration  TO  from  one  of  two  equa¬ 
tions  selected  on  the  basis  of  the  scaled  overpressure 
level. 

4.  Compute  optimum  height  of  burst  RB,  the  reciprocal  of 
the  shock  velocity  UR  and  the  peak  wind  velocity  VO, 

the  latter  two  using  the  Mach  reflection  region  equations. 

5.  Compare  the  height  of  burst  HOB  with  RB  to  determine 
reflection  region.  If  the  regular  reflection  region 
applies  then  compute  height  of  burst  coordinates  RC  and 
RD,  a  correction  factor  CF  and  the  incidence  angle  THETA. 
Recompute  the  reciprocal  of  the  shock  velocity  UR  and  the 
peak  wind  velocity  VO. 

6.  Compute  the  critical  air  displacement  values  RO  and  R00. 
Compute  the  maximum  height  of  the  dust  cloud,  HST. 

Control  is  then  returned  to  the  subprogram  DUST4. 

A . 1 . 8  Output  Subprogram 

This  subroutine  contains  the  output  statements  for 

the  off-line  listing  of  the  variables  computed  by  DUST1, 

DJST2  or  DUST3  subroutines.  The  values  of  RUNTYP,  the  type 

of  dust  computation,  and  NB,  the  number  of  bursts  are  used 

to  select  the  appropriate  output  statements.  In  addition, 

an  execution  terminated  message  is  listed  if  the  value  of 

CUTOFF  is  false  as  provided  for  in  each  of  the  dust 

subroutines . 

The  output  provides  the  spectral  dust  density 
variation  as  well  as  the  total  dust  density  as  functions 
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of  time  for  each  elevation  selected.  The  major  computed 
output  variables  are  the  observation  times,  TOBS(I),  the 
total  dust  density,  DY(T,L),  and  the  dust  densities  per 
particle  size  class,  DYC(I,N,L),  for  N=1  to  NCLASS,  L=1 
to  NY,  and  1=1  to  IMAX  in  steps  IOUT. 

A. 1.9  Uniform  Ground  Air  Temperature  Subprogram 

An  air  temperature  subroutine  was  written  to 
compute  the  air  temperature  and  dust  density  at  the  ele¬ 
vation  of  interest  for  a  site  which  is  surrounded  by  a 
uniform  ground  condition.  This  subroutine  is  entered 
through  a  CALL  ATI  statement  in  the  main  program.  Execu¬ 
tion  then  proceeds  through  three  parts:  (1)  an  initial 
computation  in  which  a  variety  of  numbers  are  generated 
for  subsequent  use,  (2)  initialization  of  dimensional 
statements,  and  (3)  the  main  calculation  loops.  The  first 
part  consists  of  the  following  steps: 

1.  Compute  initial  cell  size,  DY. 

2.  Compute  the  number  of  cells  needed,  JMAX. 

3.  Compare  JMAX  with  storage  allotment  for  J  (if  exceeded 
exit  from  subroutine) . 

4.  Compute  depth,  ZA(M) ,  in  soil  where  soil  temperature  is 
to  be  evaluated.  Set  ZM(10)  =  5. 

5.  Call  the  soil  temperature  subroutine,  ST (ALPHA,  FK, 
TEMPA)  . 

6.  Compute  air  column  motion  parameters,  RO  and  TOM. 

7.  Compute  weight  of  air  in  cell,  WAA. 

8.  Compute  depth  of  soil  eroded  per  time  step,  DZ,  after 
evaluating  TP  and  VP. 

9.  Compute  weight  of  soil  eroded  per  time  step,  W. 

10.  Compute  critical  index  for  erosion  cutoff,  IC  and 
FLOAT  (IC-1). 

The  second  part  consists  of  the  following  initialization: 

1.  The  observation  time,  T0BS(1). 

2.  The  dust  density,  DDY(l). 

3.  The  conditions  for  the  ground  level  variables,  TEY(l)  , 
YY(1>,  S(l) ,  Q(l). 
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4.  The  soil  depth  index,  M. 

5.  The  values  of  T,  TT,  PI. 

6.  The  air  temperature  variable  adjacent  to  the  ground, 

TE(1) . 

7.  The  cell  variables  at  the  shock  front,  S>JMtf(J), 

SUMWP(J),  TE(J) ,  and  S(J). 

The  third  part  consists  of  the  following  calculations: 

1.  Compute  the  value  of  constants  T1  and  G. 

Compute  maximum  time  of  interest,  TOUT. 

2.  DO  loop  on  I,  the  time  variable, 

a.  transfer  cell  temperatures  TEJ(J), 

b.  compute  FI  and  FI1, 

c.  compute  height  of  cloud,  H.  Compare  with  cloud 
height  at  end  of  positive  phase  duration.  If  H 
is  greater,  set  T=T1, 

d„  using  the  above  cloud  height  test,  select  one  of 
two  equations  with  which  to  evaluate  the  total 
elapsed  time,  T.  A  Newton-Ralphson  iteration 
procedure  is  used  and  an  accuracy  test,  based 
upon  0.0001*TM,  is  made  on  the  convergence  step, 

e.  again  using  the  above  cloud  height  test  select 
one  of  two  equations  with  which  to  evaluate  the 
observation  time,  TOBS(I), 

f.  compute  absolute  pressure,  P2,  pressure  change,  DP, 
and  a  related  quantity,  TERM2, 

g.  compute  depth  of  erosion,  Z.  Compare  Z  with 
current  temperature  depth  interval.  If  exceeded 
index  M  or  otherwise,  compute  current  value  of 
soil  temperature  TE(I), 

h.  compute  total  quantity  of  dust  in  air  column, 

WA(I). 

3.  DO  loop  on  J,  the  height  variable  (this  is  a  nested 
DO  loop) , 

a.  compute  FJ, 

b.  determine  which  cloud  region  is  applicable,  then 
proceed  to  one  of  two  groups  of  statements  to 
evaluate  quantity  of  dust  in  cell,  SUMW(J), 

c.  compare  J  with  2  to  determine  if  cell  is  adjacent 
to  soil,  then  use  one  of  two  equations  to  compute 
the  quantity  of  dust  entering  the  cell,  SUMWP(J), 

d.  compute  the  quantity  of  dust  in  air  column  above 
cell.  Evaluate  thermal  radiation  flux  incident 
on  cell,  Q(J),  and  then  compute  the  associated 
energy  deposited  in  cell, 

e.  compute  new  air  temperature  in  cell,  TE(J),  after 
evaluating  quantities  DTN  and  DTE.  Also  evaluate 
TERM1. 
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f.  compute  cell  size,  S(J),  and  cell  center  height. 
YY(J) . 

g.  Compute  dust  density  in  cell,  DD(J).  (End  of 
J  DO  loop) . 

4.  Continue  in  I  DO  loop  to: 

a.  reinitialize  T1  and  PI, 

b.  evaluate  temperature  TEY(I)  and  dust  density  DDY(I) 
at  the  height  YC  by  linear  interpolation.  This  is 
done  after  YC  is  sequentially  compared  with  the  cell 
height  YY(J).  A  quantity  DYY  is  then  computed. 

(End  of  I  DO  loop) . 

Control  is  then  returned  to  the  main  program  for  subsequent 
execution. 

The  major  computed  output  variables  are  the  obser¬ 
vation  times,  TOBS(I),  the  air  temperatures,  TEY(I),  and  the 
dust  densities,  DDY(I),  for  1=1  to  IMAX  in  steps  IOUT. 

A. 1.10  Bi-Uniform  Ground  Air  Temperature  Subprogram 

An  air  temperature  subroutine  was  written  to  com¬ 
pute  the  air  temperature  and  dust  density  at  the  elevation 
of  interest  for  a  site  which  consists  of  two  uniform-ground- 
condition  regions  in  front  of  the  site.  This  subroutine  is 
entered  through  a  CALL  AT2  statement  in  the  main  program. 
Execution  then  proceeds  through  three  parts:  (1)  an  initial 
computation  in  which  a  variety  of  numbers  are  generated  for 
subsequent  use,  (2)  initializing  of  dimensioned  statements, 
and  (3)  the  main  calculation  loops.  The  first  part  consists 
of  the  following  steps: 

1.  Compute  initial  cell  size,  DY. 

2.  Compute  the  number  of  cells  needed,  JM4X. 

3.  Compare  JMAX  with  storage  allotment  for  J  (if  exceeded, 
exit  from  subroutine). 

4.  Compute  depth,  ZA(M),  in  soil  where  soil  temperature  is 
to  be  evaluated.  Set  ZM(10)  =  0. 

5.  Call  the  soil  temperature  subroutine  ST(ALPHA,  FK,  TEMPA). 

6.  Compute  air  column  motion  parameters,  RO  and  TOM. 

7.  Compute  weight  of  air  in  cell,  WAA. 
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8.  Compute  depth  of  soil  eroded  per  time  step,  DZ,  after 
evaluating  TP  and  VP. 

9.  Compute  weights  of  soil  eroded  per  time  step  for  each 
zone,  W1  and  W2. 

10.  Compute  critical  index  for  range  cutoff,  ISTAR. 

The  second  part  consists  of  the  following  initializations: 

1.  Air  temperature  at  ground  level,  TEY(l). 

2.  Observation  time,  T0BS(1). 

3.  Conditions  for  the  ground  level  variables  DDY(l),  YY(1), 

S (1) ,  and  Q(l) . 

4.  Logical  value,  QQ. 

5.  The  values  of  Tl,  V,  and  PR1. 

6.  The  air  temperature  variable  adjacent  to  the  ground, 

TE(1) . 

7.  The  cell  variables  at  the  shock  front,  SUMW(J) ,  SUMWP(J), 
TE(J) ,  and  S(J). 

The  third  part  consists  of  the  following  calculations: 

1.  Compute  the  number  of  integration  paths  to  be  taken, 

NMAX.  Compare  this  with  storage  allotment  for  N.  If 
exceeded,  set  NI4&X  equal  to  maximum  permissible  value. 

2.  DO  loop  on  N,  the  observation  time  variable, 

a.  evaluate  in  which  zone  integration  starts.  If 
integration  starts  in  zone  1,  evaluate  critical 
indices  IS,  IM,  IC  and  proceed.  If  integration 
starts  in  zone  2  and  if  this  is  not  the  firt  time 
this  occurs,  evaluate  critical  index  IS  and  update 
M,  Tl,  V,  PR1,  TE(1) ,  SUMW(J),  SUMWP(J),  TE(J), 
and  S(J).  If  integration  starts  in  zone  2  for 
the  first  time,  change  logical  value  QQ  and 
initialize  M,  Tl,  V,  PR1,  TE(1),  SUMW(J) ,  SUMWP(J), 
TE(J),  and  S(J).  If  integration  starts  in  zone  2, 
evaluate  IC  and  IM, 

b.  evaluate  FC. 

3.  DO  loop  on  I,  the  transit  time  variable  (this  is  a 
nested  DO  loop) . 

a.  transfer  cell  temperatures  TEJ(J), 

b.  compute  FI  and  HI, 

c.  compute  velocity  correction  VC,  new  time  T2,  and 
velocity  V. 

d.  compute  observation  time,  TOB, 

e.  compute  absolute  pressure,  PRZ,  pressure  change, 

DP,  and  a  related  quantity,  TERM2. 
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f.  compute  depth  of  erosion  Z.  Compare  Z  cuLj.er.t 

temperature  depth  interval.  If  exceeded,  index  M 
and/or  otherwise  compute  current  value  of  soil 
temperature  TE(I), 

g.  compute  total  quantity  of  dust  in  air  column,  WA(I), 
from  one  of  two  formulas  depending  upon  starting 
zone , 

h.  if  W.A(l)  is  very  small,  set  Q1  "  0:  otherwise  set 
Q1  =  Q0. 

4.  DO  loop  on  J,  the  height  variable  (this  is  a  double 
nested  DO  loop) . 

a.  compute  FJ, 

b.  determine  which  cloud  region  is  applicable,  then 
proceed  to  one  of  two  groups  of  statements  to 
evaluate  quantity  of  dust  in  cell,  SUMW(J)t 

c.  compare  J  with  2  to  determine  if  cell  is  adjacent 
to  soil,  then  use  one  of  two  equations  to  compute 
the  quantity  of  dust  entering  the  cell,  STJMirfP(J)f 

d.  compute  the  quantity  of  dust  in  air  column  above 
cell.  Evaluate  thermal  radiation  flux  incident 
on  cell,  Q(J),  and  then  compute  the  associated 
energy  deposited  in  cell, 

e.  Compute  new  air  temperature  in  cell,  TE(J),  after 
evaluating  quantities  DTN  and  DrE.  Also  estimate 
TERM1, 

f.  compute  cell  size,  S(J),  and  cell  center  height, 
YY(J), 

g.  compute  dust  density  in  cell,  DD(J).  (End  of 
J  DO  loop) . 

5.  Continue  in  I  DO  loop  to: 

a.  reinitialize  T1  and  PI, 

b.  evaluate  temperature  TEY(l)  and  dust  density  DDY(I) 
at  the  height  YC  by  linear  interpolation.  This 

is  done  after  a  YC  is  sequentially  compared  with 
the  cell  height  YY(J).  A  quantity  DYY  is  then 
computed.  (End  of  I  DO  loop). 

Control  is  then  returned  to  the  main  program  for  subsequent 

execution. 


The  major  computed  output  variables  are  the  obser¬ 
vation  times,  TOBS(I),  the  air  temperatures,  TEY(I),  and  the 
dust  densities,  DDY(I),  for  1=1  to  IMAX  in  steps  IOUT. 

A. 1.11  SOTEMP  Subprogram 

This  subroutine  was  written  to  compute  the  tem¬ 
perature  distribution  in  the  soil  at  the  time  of  shock 
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arrival  resulting  from  the  thermal  radiation  flux  of  the 
fireball  of  a  nuclear  burst.  Reradiation  at  the  surface 
is  also  included  in  the  computations.  This  subroutine  is 
entered  through  a  CALL  SI (ALPHA,  FK,  TEMPA)  statement  in 
either  of  the  air  temperature  subroutines.  Execution  then 
proceeds  through  three  basic  parts:  (1)  an  initial  compu¬ 
tation  in  which  a  variety  of  numbers  are  generated  for 
subsequent  use,  (2)  computation  of  the  temperature  and 
and  reradiation  at  the  soil  surface,  and  (3)  computation 
of  the  temperature  distribution  and  gradient  below  the 
surface.  The  first  part  consists  of  the  following  steps: 

1.  Compute  the  values  of  TEMPAS,  TEMPA4 ,  GAMMA,  C0NS1, 

CONS 2. 

2.  Compute  the  time  step  size,  DELT. 

3.  Compute  the  values  of  CONS 3  and  C0NS4. 

4.  Compute  the  times,  T(M). 

5.  Compute  the  incident  thermal  radiation  flux  rate,  QI(M). 
The  second  part  consists  of  the  following  surface -condition 
steps: 

1.  DO  loop  on  M,  the  time  variable:  at  each  time  compute 
reradiation  flux  rate,  QR(M),  the  net  radiation  flux 
rate,  G(M),  and  the  piece-wise  linear  G(M)  curve  param¬ 
eters,  A  (Ml)  and  B(M1). 

2.  Inner  DO  loop  on  I,  the  dummy  time  integration  variable, 

a.  sum  the  terms  corresponding  to  time  steps  prior  to 
M-l  that  contribute  to  the  temperature  at  M;  com¬ 
pute  SUM, 

b.  compute  by  iteration  the  temperature  at  the  surface 
corresponding  to  M;  compute  TEMP. 

The  third  part  consists  of  the  following  temperature  distri¬ 
bution  steps: 

1.  Set  TS(1)  equal  to  the  last  computed  value  of  TEMP  from 
the  second  part. 

2.  DO  loop  on  N,  the  soil  depth  variable:  compute  the 
values  of  Cl,  C2,  C3,  C4,  El,  and  EF1. 

3.  Inner  DD  loop  on  I,  the  dummy  time  integration  variable, 

a.  compute  the  values  of  C5,  C6,  C7,  II,  E2,  and  E3, 

b.  compute  the  temperatures,  TS(N). 
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Control  Is  then  returned  to  the  calling  subroutine  for 
subsequent  execution. 

A. 1.12  FIRBAL  Subprogram 

This  is  a  function  subprogram  named  FUNCTION 
PE(TM4X,T)  where  the  argument  T  provides  the  current  time 
at  which  a  value  is  to  be  assigned  to  PE.  The  scaled 
fireball  thermal  power  is  compiled  directly  into  the  sub¬ 
program  by  a  table  of  values  of  scaled  time  and  scaled 
power  listed  in  data  statement  format.  The  scaled  time 
corresponding  to  T  is  computed,  and  then  the  value  of  PE 
is  evaluated  by  a  standard  linear  interpolation  formula. 

This  subprogram  is  entered  when  its  name  is 
referenced  in  an  arithmetic  statement  of  the  calling  sub¬ 
program,  i.e.,  SOTEIMP.  Control  and  the  value  of  PE  are 
then  returned  to  said  arithmetic  statement. 

A. 2  Sample  Problem 

Sample  computer  prediction  results  were  obtained 
for  the  dust  density  and  air-temperature  environment  at  a 
sample  site.  The  computer  job  consisted  of  three  data  sets 
to  illustrate  the  computer  program  operation  and  provide 

•  uniform  ground,  single-burst  (D'JSTl)  dust 
computations , 

•  nonuniform  ground,  single-burst  (DUST2) 
dust  computations, 

•  uniform  ground,  single-burst  (ATI)  air- 
temperature  computations. 

The  nuclear  attack  condition  assumed  was  one  20-MT 
weapon,  detonated  at  optimum  height  of  burst,  with  a  site 
peak  overpressure  of  10  psig. 

The  sample  site  was  located  on  a  quadrangle  map; 
an  overlay  map  was  made  of  the  area  in  and  around  the  site 
location.  The  overlay  map  features  distinctive  character¬ 
istics  of  the  locale  that  would  be  influential  to  the  dust 


A-31 


and  air-temperature  environments  at  the  site.  Locations  of 
the  site  installations,  radar,  missiles,  and  so  forth,  were 
also  drawn  on  the  overlay  map.  The  overlay  map  for  this 
sample  problem  is  illustrated  in  Figure  A.l. 

The  numbers  shown  on  the  overlay  map  correspond 
to  the  condition  of  the  soil  surface  in  that  area.  One 
might  expect  that  different  land  usages  would  have  varying 
effects  on  the  soil  particles  that  are  able  to  rise.  Each 
of  IITRI's  dust  models  include  weight  fraction  terms  for 
the  erodible  soils.  Under  varying  land  conditions,  all  of 
the  soil  may  not  be  erodible  and  assumptions  have  to  be 
made  as  to  the  fraction  of  the  soil  that  is  erodible. 

Table  A. 4  is  a  legend  of  the  land  usage  factors  used  for 
various  land  usages.  Table  A. 5  describes  the  particle  size 
distribution  and  atmospheric  variables  of  the  sample  site. 

The  soil  class  numbers  correspond  to  the  soil  classifications 
as  tabulated  on  Table  A. 6. 

A. 2.1  Nonuniform  Ground  Dust  Model  Variables 

Site  input  variables  were  determined  for  the  dust 
computations  for  the  sample  site  in  the  following  manner. 

It  was  assumed  that  the  attack  azimuth  is  the  same  as  the 
primary  orientation  of  the  radar.  A  line  was  drawn  parallel 
to  the  radar  orientation  and  through  the  location  of  the 
intake  and  exhaust  ports  (+) .  The  land  over  which  this  line 
crosses  within  the  3000-ft  radius  is  the  soil  surface  that 
will  contribute  to  the  dust  loading  at  the  intake  and  exhaust 
ports.  The  land  along  this  line  was  divided  into  land  usage 
classifications.  In  front  of  the  intake  and  exhaust  ports, 
the  land  usage  is: 

0  to  200  ft  --  Nonerodible  land  (1),  zero 

soil  particle  contribution  to 
the  dust  cloud. 
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Table  A. 4  LEGEND  OF  LAND  USAGE  FACTORS 


r  ■  •  : — — — r— —  : — — — — — . .  — -  ■  i 

Land 

Classification 

Land  Description 

Land  Usage 
Factors 

1 

Launch  areas,  concrete  roads, 
radiation  perimeter,  barren  rock, 
200-ft  diameter  around  intake  j 

and  exhaust  ports,  etc. 

0  erodible 
soil 

2 

| 

Thick  grassland 

0.20  of  soil 
classes  1 
and  2 

3 

Scrubland  and  woodland  areas. 

0.50  of  all 
soil  classes 

4 

Farmland,  cultivated  areas, 
orchards,  beaches,  etc. 

1.0  of  all 
soil  classes 

5 

Railroads,  gravel  pits,  etc. 

1.0  of  all 
soil  classes 

6 

Residential  and  industrial 
areas . 

0.10  of  all 
soil  classes 

7 

Water . 

0  erodible 
soil 

8 

Undeveloped  land  as  not  clas¬ 
sified  above. 

0.75  of  all 
soil  classes 

i  . 
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Table  A. 5  SAMPLE  SITE  VARIABLES 


County:  Typical 

Temperature:  ave.  max.  74°F 

yearly  ave.  50°F 
ave.  min.  25 °F 

Site  Elevation:  600  feet 

Ave.  Barometric  Pressure:  14.4  psia 

Mechanical  Soil  Analysis: 


Soil  Class 
Particle  Size 


(mm) 


'  8 

7 

6 

5 

D 

3 

2 

1 

5.0 

m 

O 

• 

o 

0.02 

.005 

0 

9 

53 

□ 

B 

Percent 


Table  A. 6  SOIL  CLASSIFICATION 


200  to  3000  ft  —  Undeveloped  land  (8),  0.75 

of  the  soil  particles  of  all 
class  sizes  will  be  contributed 
to  the  dust  cloud. 

To  the  rear  of  the  intake  and  exhaust  ports  the  land  usage  is 

0  to  200  ft  --  Nonerodible  land  (1),  zero 

soil  particle  contribution  to 
the  dust  cloud. 

200  to  1250  ft  --  Undeveloped  land  (8),  0.75 

of  soil  particles  of  all  class 
sizes  will  be  contributed  to 
the  dust  cloud. 

1250  to  1350  ft  --  Nonerodible  (1),  toll  road, 

zero  soil  particle  contribution 
to  dust  cloud. 

1350  to  1800  ft  --  Undeveloped  land  (8),  0.75 

of  the  soil  particles  of  ail 
class  sizes  will  be  contributed 
to  the  dust  cloud. 

1800  to  2300  ft  —  Woodland  (3),  0.50  of  the  soil 

particles  of  all  class  sizes 
will  be  contributed  to  the 
dust  cloud. 

2300  to  2700  ft  —  Water  (7),  zero  soil  particle 

contribution  to  the  dust  cloud. 

2700  to  3000  ft  —  Railroad  (5)  1.00  of  soil 

class  size  (8)  will  be  con¬ 
tributed  to  the  dust  cloud. 

The  particle  size  distribution  for  the  site  is  multiplied  by 
these  land  usage  factors.  The  product  is  an  adjusted  weight 
fraction  of  the  soil  available  for  erosion  in  each  range  area 

A. 2. 2  Uniform  Ground  Hast  Model  Variables 

Site  input  variables  are  determined  for  the  IITRI 
uniform  soil  condition  dust  model  in  the  following  manner. 
First,  no  attack  azimuth  need  be  chosen  as  the  soil  con¬ 
ditions  over  the  entire  3000-ft  radius  will  be  averaged; 
consequently,  soil  availability  is  the  same  for  all  attack 
azimuths.  The  ground  conditions  are  averaged  by  determining 
the  percent  of  land  area  that  is  within  the  3000-ft  radius, 
that  is,  in  each  of  the  land  usage  classifications  as 
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defined  in  Table  A. 4.  For  the  site  example  chosen  (Fig¬ 
ure  A.l),  the  ratio  of  a  particular  land  usage  area  to  the 
total  area  within  the  30D0-ft  radius  is  shown  in  Table  A. 7. 
For  example,  approximately  0.04  (4  percent)  of  the  land  is 
woodland  area  (land  usage  classification  3),  The  land  usage 
factor  for  woodland  is  0.50,  from  Table  A. 4.  This  factor 
multiplied  by  the  area  ratio  gives  0.04  x  0.50  *  0.020,  the 
soil  multiplier.  The  particle  size  distribution  for  each 
particle  size  is  then  multiplied  by  the  soil  multiplier  to 
obtain  a  weight  fraction  for  each  particle  size  within  the 
land  classification.  The  soil  distribution  at  the  sample 
site  for  particle  size  4  is  0.53.  Multiplying  the  soil 
distribution  by  the  soil  multiplier  gives  0.53  x  0.02  * 
0.0106.  This  number  represents  the  weight  fraction  of 
class  4  particles  contributed  by  the  woodland  areas.  The 
weight  fraction  for  each  particle  size  is  totaled  over  all 
land  classifications.  The  numbers  shown  at  the  bottom  of 
Table  A. 7  are  the  adjusted  weight  fractions  for  each  par¬ 
ticle  size  for  the  sample  site. 

A. 2. 3  Sample  Site  Data  Input 

The  data  input  is  provided  according  to  the  dis¬ 
cussion  in  Section  A. 1.1.  There  are  three  data  sets  cor¬ 
responding  to  the  three  types  of  computations  performed 
for  the  sample  sice.  Table  A. 8  lists  the  values  of  the 
variables  required  as  input  by  the  computer  program.  The 
FORTRAN  names  of  all  the  variables  are  given  in  parentheses 
where  they  first  appear  in  the  table. 


Table  A.  7  AVERAGE  SOIL  PROPERTIES  FOR  SAMPLE  SITE 


Table  A. 8  SAMPLE  SITE  DATA  INPUT 
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Table  A. 8  (continued) 
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Appendix  C 

WEAPONS  EFFECTS  VARIABLES 

The  dust  and  air-temperature  prediction  models 
depend  upon  certain  weapons  effects  variables.  This  appendix 
briefly  summarizes  the  variables  of  interest  and  the  empirical 
equations  which  were  developed  for  computer  application. 

The  attack  situation  is  defined  in  terms  cf  the  weapon  yield, 
the  overpressure,  and  the  height  of  burst  or  its  equivalent 
(i.e.,  scaled  height  of  burst,  optimum  height  of  burst,  or 
overhead  burst  are  the  equivalents  considered) .  The  ambient 
air  pressure  and  temperature  are  also  used  to  define  the  wea¬ 
pons  effects  environment .  The  general  range  of  overpressures 
consic'ared  herein  is  from  approximately  5  psi  to  50  psi. 

The  accuracy  level  to  which  the  pertinent  variables  are 
determined,  by  the  simplified  equations,  is  approximately 
10  percent. 

The  weapons  effects  variables  that  must  be  evaluated 

are: 

•  peak  air  velocity,  V  ,  fps, 

•  positive  phase  duration  of  wind,  tQ,  sec, 

•  exponential  constant  for  air  velocity,  k, 

•  shock  velocity,  U,  fps, 

•  shock  air  density,  pg,  lb/ft  , 

•  shock  air  temperature,  T  ,  °F, 

3 

•  average  air  density,  p  ,  Ib/ft  , 

a 

•  positive  phase  duration  of  overpressure,  t'Q,  sec, 

•  overpressure  waveform  parameter  A,  a,  and  3, 

•  time  of  arrival  at  site  t  ,  sec, 

cl 
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•  time  of  maximum  thermal  power,  t  ,  sec , 

2 

•  maximum  thermal  power,  P_  .  Btu/ft  sec, 

r  1  max 

•  atmospheric  attenuation,  t,  and 

•  the  angle  of  incidence,  0. 


The  weapons  effects  data  were,  for  the  most  part,  obtained 

C  1 

from  Gladstone.  * 


The  scaled  height -of -burst  for  overpressure  were 
fitted  with  two  straight  line  segments:  one  line  connecting 
the  overhead  burst  point  with  the  optimum  height-of -burst 
point,  and  the  second  line  connecting  the  optimum  height-of - 
burst  point  with  the  surface  burst  point .  These  two  straight 
line  segments  are  defined  by  four  scaled  quantities:  Hg(j  is 
the  scaled  overhead  burst  height,  is  the  scaled  optimum 

height  of  burst,  R  is  the  scaled  range  for  optimum  height 
of  burst,  and  R  is  the  scaled  range  for  a  surface  burst. 

Set 

The  scaled  overpressure  p  is  defined  as 

s  & 


p 


ss 


(C.l) 


where 

p  =  actual  overpressure, 

pgo  =  standard  ambient  pressure,  14.7  psia,  and 
PQ  =  actual  ambient  pressure. 


This  approximate  scaling  law  is  used  because  it  is  expected 

that  P  will  deviate  only  slightly  from  P  .  The  scaled 
o  so 

height  of  burst  and  ranges  are  defined  as 


- 


H 


-s  WT73 


so 


1/3 


(C.2) 


where 


H  *=  the  scaled  height  or  range,  ft/MT1^, 
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H  =  the  corresponding  actual  height  or  range, 
ft,  and 

W  =  weapon  yield,  MT . 

The  four  scaled  quantities  defining  the  height -of -burst 
approximations  are: 


(C.3) 

(C.4) 

(C.5) 

(C.6) 


Another  scaled  height-of -burst  factor  will  be  needed  for  the 
time-of -arrival  calculation.  This  scaled  quantity  is 


Hse=7<Hsd  +  Hsb>‘  (C*7) 

A  number  of  the  weapons  effects  variables  depend 
only  upon  the  scaled  overpressure  psg.  These  relate  primarily 
to  the  wind  and  overpressure  waveforms.  The  exponential  coef¬ 
ficient  k  used  in  the  velocity  waveform  has  been  approximated 

as 

k  =  0.877  pss  -  0.3,  (C,8) 


with  a  minimum  value  of  k-1  for  the  lower  overpressure  range. 
The  overpressure  waveform  parameters  were  obtained  from 
Erode . ^ ^  A  double  exponential  waveform  is  used  for  the 


overpressure  range  of  interest',  that  is 


ps(t>  =  p. 


-  e‘a(t/to^  +  B  e‘P*t/fto)j 


where 


(C.9) 


t  ■  time, 

t^  =  positive  phase  duration,  and 
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B,  a,  and  3  are  the  waveform  parameters.  The  following 
empirical  equations  were  developed  for  these  parameters: 


B  =  0.092  ps°-328  , 

(C.10) 

/V  n  OQQ  0.318  J 

a  =  0.398  p  ,  and 

S  S 

(C.ll) 

fi  —  °  1  _  0*262 

P  ■  -'-13  Pss 

(C. 12) 

The  positive  phase  duration  of  the  wind  and  the 
overpressure  are  both  relatively  independent  of  the  height 
of  burst  and  can  therefore  be  formulated  in  terms  of  the 


overpressure  alone.  The  time  variables  scale  in  the  follow 


ing  manner: 

t 


(C.  13) 


where 

t  =  actual  time  variable, 
t  =  scaled  time  variable, 

Tos  -  standard  ambient  temperature ,  70°F,  and 
Tq  =  actual  ambient  temperature. 

Temperature  must  be  specified  by  its  absolute  value. 

The  following  empirical  equations  were  developed: 


tos  =  3.59  -  0.107  pss  pss  ^  30  psi  (C.14) 
or  tQS  =  4.125  -  °.°285  pss  pss  2  30  psi  (C.15) 

where  tos  is  the  positive  phase  duration  of  the  velocity. 

t'  “  2-5  -  0.024  poc  +  1.5  e'0,16  pss  (C.16) 

OS  ss 

where  t*  is  the  positive  phase  duration  of  the  overpressure. 
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The  value  of  some  of  the  weapons  effects  variables 
depend  upon  both  the  strength  of  the  incident  wave  and  the 
nature  of  the  shock  reflection  process.  For  the  range  of 
overpressures  of  interest  and  the  accuracy  requirements,  the 
influence  of  the  shock  reflection  process  can  be  eliminated 
for  some  of  the  variables.  Thus,  the  temperature  and  density 
of  the  air  at  the  shock  front  can  be  determined  from  the  con¬ 
ventional  Rankine-Hugoniot  equation.  The  average  air  density 
has  been  determined  as  a  simple  arithmetic  mean  between  the 
values  at  the  shock  front  and  the  ambient  density. 

The  apparent  shock  velocity  (i.e.,  the  velocity  at 
which  initial  pressure  disturbances  sweep  over  the  ground 
surface)  and  the  horizontal  component  of  the  air  velocity  at 
the  shock  front  are  strongly  dependent  upon  the  shock  config¬ 
uration.  The  shock  configuration  can  be  divided  into  two 
regions:  (1)  the  Mach  reflection  region,  and  (2)  the  regular 
reflection  region.  These  shock  reflection  configurations  are 
illustrated  in  Figure  C.l. 

The  Mach  reflection  process  will  exist  for  all 
heights  of  burst  below  the  optimum;  that  is,  the  optimum 
height  of  burst  is  the  approximate  transition  point  between 
the  two  reflection  regions.  The  shock  velocity  and  the  peak 
air  velocity  are  both  given  by  classical  Rankine-Hugoniot 
relations.  The  time  of  arrival  of  the  shock  wave  does, 
however,  vary  with  the  hieght  of  burst  within  this  region. 

The  following  empirical  equation  was  developed  for  the 
scaled  time  of  arrival. 

tas  -  79.0  Ps;0,88  (0.4  +  0.6  Hs)  (0.17) 

The  reflection  process  in  the  regular  reflection 
region  has  been  idealized  as  simple  linear  process  such  that 
the  angle  of  incidence  a  is  equal  to  the  angle  of  shock 
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(a)  Mach  Reflection  Region 


(b)  Regular  Reflection  Region 

Figure  C.l  Shock  Configuration  at  Ground  Surface 
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reflection  and  that  the  strength  of  the  two  waves  are  equal. 
The  pressure  is  known  in  the  reflection  region  (State  sp. 
These  assumptions  permit  us  to  compute  the  pressure  behind 
the  incident  wave  (State  Sp.  Thus,  the  shock  velocity  Uj 
and  the  particle  velocity  associated  with  the  incident 
wave  are  known.  The  following  approximate  relations  were 
thus  obtained: 


U{  ■=  c0  Vo. 572  +  .428  ? 


(C.18) 


and 


0.945(g-l)cp 
*Y4  +  3  ? 


where 

§  *  overall  pressure  ratio  =  ps/PQ,  and 

cQ  =  ambient  sound  velocity. 

The  apparent  shock  velocity  U  is  given  as 

ui 


(C.19) 


(C.20) 


where  a  is  the  angle  of  incidence.  Also,  the  horizontal 
component  of  the  air  velocity  is 

V-l  »  ?  V{  .  (C.21) 


The  use  of  the  above  approximate  regular  reflection  equations 
results  in  values  of  both  shock  velocity  and  particle  veloc¬ 
ity,  at  the  optimum  height  of  burst  point,  which  differ 
slightly  from  the  values  obtained  from  the  corresponding 
Mach  reflection  equations.  This  discrepancy  is  of  the  order 
of  two  to  four  percent  and  is  a  function  of  the  strength  of 


p 
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the  shock  wave.  In  order  to  eliminate  this  discrepancy 
the  following  equations  were  developed  and  used. 


U  =  Ur  •  f/sin  a 


(C. 22) 


Vl  =  Vr  sin  a/f  (C. 23) 

where 

f  =  sin  (arctan  (Rsc/Hsb>),  (C.24) 

U  -  Kankine-Hugoniot  shock  velocity,  and 
Vr  =  Rankine-Hugoniot  particle  velocity. 


The  scaled  time  of  arrival  in  the  regular  reflection  region 
is  calculated  from  the  following  empirical  equations: 


and 


as 


<t2-tl>(Hse-Hs> 

^  +  — W  ’  "s 


se  ’ 


30  Pss-0672  . 


•7Q  „  -0.88 

/9  Pss 


(0.25) 

(C. 26) 
(0. 27) 


30 


-0.672 


(C.  28) 
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